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ABSTRACT 


As part of the SEASAT program of NASA, two sets of 
analysis programs were developed for ECON, Inc. under 
contract NASW-21358 . One set of programs produce 63 x 63 
horizontal mesh analyses on a polar stereographic grid. 

The other set produces 187 x 187 third mesh analyses. 

The parameters analyzed include sea surface temperature, 
sea level pressure and twelve levels of upper air temperature; 
height and wind analyses. Both sets use operational data 
provided by Fleet Numerical Weather Central. The analysis 
output is used to initialize the primitive equation fore- 
cast models also included as part of this contract. 
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SECTION I. SCALAR ANALYSIS USING THE PATTERN -CON SERVING 
TECHNIQUE l 

A. Introduction 

In meteorological analysis, one piece of information 
that should not be ignored is the most recent past analysis, 
or, if available, a good forecast valid at the current ana- 
lysis time. A man doing a hand analysis usually uses such 
information. In particular, the analyst needs an estimate 
of the positions of highs and lows (curvature) and of the 
areas of strong and weak gradients. In referring to the 
past analysis or forecast, the man usually is looking not 
so much for absolute magnitudes as for the shape of the 
field. When he draws his new analysis, he will first attempt 
to fit the shapes in the' past field to the new data. If 
conflict occurs, the new data takes precedence unless the 
analyst suspects that the data is in error. In regions;; where 
data is routinely sparse, many conflicts need to be resolved 
because of the accumulation of errors. This results from : 
the cycle of a poor analysis initializing a forecast which 
is consequently poor which is used as a deficient first 
guess for the next analysis/prediction cycle. 

The best objective analysis scheme is probably one 
that follows the same rules that a man follows. If such 
objective techniques can be worked out, the machine may do 
the job in a more consistent manner than a man is able to do. 
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The basic goal of this analysis is to fit the follow- 
ing information to varying degrees: the new data; the most 

recent past analysis or forecast value (the first guess); 
the gradients of the first guess in eight directions from 
each grid point; and the Laplacian of the first guess. 

The degree of fit desired for each piece of information is 
Specified by an array of weights. These variables and 
weights are named in Table .1-1. 

The desired fit is realized by minimizing the sum of 
the deviations of the various characteristics of the ana- 
lysis from their counterparts in the first guess. The 
minimization is accomplished with an elementary application 
of the calculus of variations. 

Information is spread through space by the gradient 
and Laplacian terms. In a surface analysis, there are 
sometimes natural obstacles (mountain ridges, coastlines, 
etc.) beyond which an analyst would not allow a new obser- 
vation to influence the analysis. This kind of constraint 
can be simulated in the objective analysis by reducing the 
weights of the gradients and Laplacian along the demarcation 
zone. 

The decision on the magnitudes of the various weights 
is less arbitrary if we view each weight as the inverse of 
the variance associated with the parameter it multiplies. 
This viewpoint is also useful in re-evaluating the weights 
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of the data. 


An analysis cycle consists of three basic steps: 

1 . Assemble the data at grid points . 

2. Solve the minimization equation. 

3. Re-evaluate the weight of each report. 

In order to adequately evaluate the weight of each report, 
at least two cycles are required. It is desirable to in- 
clude one additional cycle to allow initially suppressed 
data to enter the analysis with a high weight if supporting 
data some distance away causes the analysis to conform 
more closely to the report after the second cycle. The 
basic steps are detailed individually in the following 
sections . 



B. 


Assembly 


We shall refer to the guess field as P^ with weight 

A. . . On the first cycle, it is the first guess, and A. . 

has a low and probably uniform value. On subsequent cycles, 

P, .is the result of the previous cycle, but A. . reverts 
i/3 i/3 

back to its original value. During assembly, A. . will be 

i/3 

modified to incorporate the data weights, as re-evaluated 
-at the end of the previous cycle. 

Each report is moved to the nearest grid point along 
the gradient of the guess field; that is, the guess field 
is interpolated to the report location and the difference 
between the interpolated value and the nearest grid point 
value is added to produce a pseudo reported value. A weighted 
average of the guess value and all reports affecting the grid 
point is done. 


. 

1/3 


A i,j P i,j + A 1 P 1 + A 2 P 2 + 


■ + A P 
n n 


A. . + A n + A 0 + 
1 2 


+ A 


n 


L . = A. . + A. + A~ + ... + A , 

1/ j i/3 1 2 n 


The superscript ('v) designates assembled values. A^ is 
the weight of the n report and P n is the n report. 

In the following sections, the superscript (M will not be 


used. The assembled values will be referred to as P. . » 

i/3 
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c. 


Minimizing the Deviations 


TABLE I- 1 : PCT SCALAR CONSTRAINTS 






Constraint 


Weight 

Variable being analyzed (assembled value) 

A ir j 

y axis gradient = 

P i,j4-1 " P i,j 

B i t j 

x axis gradient = 

P i+l,j P i,j 

C. . 
i,3 

x~l,y+l gradient 

= p . - p 

l-i, 3+1 1,3 

E . . 
i,3 

x+l,y+l gradient 

= P . , - P 

i+l, 3+1 i, j 

F if j 

Laplacian = P^ + ^ 

. + P. . . + p. . . + P. , ' ■ - 4p . , 

3 1-1,3 1,3+1 1,3-1 i,3 

°i , j 


The first guess shapes y, v, a, 3 and L and their res- 
pective weights B, C, E, F and D have a constant value during 
the entire analysis . Within limits specified by the weights, 
we require the final analysis to have similar values of the 
constraints as the first guess field. 
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To effect this matching, we shall minimise the following 
integral: 


(a) 

(b) 

(c) 

<d) 

(e) 

(f) 

(g> 

(h) 

(i) 

(j) 
00 
( 1 ) 

(m) 

(n) 


I = // [ 

. (P* . - P. .) 2 + 

1,3 1,3 1,3 


B i,j - n.i - "i.f + 

3 2 > 


B . . , (P* . - p* . - u 

h]-l l ( ] 1,3-1 l, j-1 

C. . (P* , . -P* . - V. . 7 

1,3 1+1,3 1,3 1,3 


) 2 + 


Z. , . (P* . - p* - V 

1-1,3 1 1,3 i-l, j l 

E. . (P* . - Pf , - a, . ) 2 + 

1,3 i-l, 3+1 1,3 i,3 

E i+l,j-l (P i,j ‘ ?|+l,j-l “ a i+l,j-l 


- V. . ) 2 + 

-1,3 


) 2 + 


) 2 + 


[1.1] 


F i,j ^+1,3+1.^!,^“ 0 ifj ) 2 - 
F i-1, j~l (P i,j ~ P i-l, j-1 ~ 6 i-l. -j-1 

D i,j .<*!«, j + + n,i+ 1 + n,i-i - *n , 3 - Vif 2 * 

Vm < p r,j + n-2>) - 4p |.y - ** 

+ p i + i,j-i - 4p - 


D 


D 


' i+1 , j 


i+i , j <n«.j - p !.j + p i+i,i+i 

<p i+i ( j-x + n.i * n,s-i 

D i,j+1 (P i+l,j+l ‘ 1 1-3 i j + 1 + P l f j+2 + P i*,j - 4p J, j + l 
Jdxdy 


- L. 


- 4P* • . - I 
1,3-1 


In the above, the starred quantities are the analysis 

values we are seeking. Each term is a departure from the 

desired matching of differential properties. Extra terms 

have been added to account for the effect of a changing P* . 

1 / J 

on the differential properties computed at surrounding points. 
Their effect is to more closely couple neighboring grid 
points. See Figure 1-1 for a depiction of the minimization 
stencil as it relates to the terms of equation [I . 1] . 

To minimize the integral, we simply take the first variation 
with respect to Pt . , and set it to zero (see equation [1.2]). 

i , j 

The solution of the resulting equation will be the Pt ^ that 
will cause the integral to be minimized. The fact that each 
term is squared insures a minimum as opposed to a maximum 
value. 
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© 


LEGEND : ( ) = constraint from equation [I. 1] 

O = difference 
>■ ~ gradient 

: laplacian at grid point 



@ = grid points 


FIGURE I- 1 : SCALAR MINIMIZATION STENCIL 
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W* = // 1 2R i , i £P 1,J - p i,:> 

~ 2 B i, j (P i, j+1 " P i,j V i,j ) 

i,3~l '1,3 i, 3-1 H i, 3’".1 / 

- 2 C. . (P* 5; . - P* . - v. .) 

if] i+l , 3 1,3 1,3 

+ 2 C. , . (P* . - P* v • " v. . ) 

1-1,3 1,3 1-1,3 1-1,3 

- 2 e m "i-i, .wi - ‘ij - “i (j > 

+ 2 E i+1 , j-1 * P i,j ” P i+1, j-l ” a i+l , j-1^ 

■ a .v; 

- 2 F. . (P* , - P* , - g. .) 

1,3 l+l, 3+1 1,3 1,3 

+ 2 F. .. . . (P* . - P* . , . . , - 0. , . . ) 

l-l, 3-1 i,3 i-l, 3-1 l-l, 3-1 

- 8 D. . (P* . + P* . . + P* x + P* . ' - 4P* . - L. . ) 

1,3 1+1,3- 1-1,3 1,3+1 1,3-1 1,3 1,3' 

+ 2 D. , . (P* . + P* „ , + >* . 

1-1,3 1,3 1-2,3 1-1,3+-' 


[ I .2] 


+ 2 D. 


. (P* 


+ P* . + P* 


+ 2 D i,i-x "ta.l-X + ’l-X.J-X 


1,3 


+ 2 D, 


(P* 


+ P* 


+ P* ' . - 

1 - 1 , 3-1 

- 4P* . . -L. , .) 

1-1,3 1-1,3 

+ p* 

i+1,3-1 

~ 4P* ,, . -L . , . ) 

1+1,3 1+1,3 

. + P* . 

4P* . , -L. . . ) 

1,3-1 1,3-1 

3 1,3-2 

. ' + P* . 

- 4P* . , ~L . . n 

,3+2 1,3 

1,3+1 1,3+1 


] dxdy s f fc 0 


• • 5 i ' 

The terms in can be grouped into three 
categories: 

1. Those involving P*. ... . . V 

■ !■;- : 1 , D 

2. Those involving P* at surrounding points 

3. Those not involving P*. 
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[A. . + B , . + B . . , + C . , + C . .. . + E . . 

1/3 i/3 i,3”l 1/3 1-1/3 1/3 

S. .P* • I + E . , . , +' F. . + F. . ' . + 16 D. .■ + D. - 

1/3 i/3 i x+l/j~l i/j i — 1/3“ 1 i/3 i“l/, 

+ D . , , . + D . . . + D . P* . 

1+1/3 i/3-l i/3+l i,3 


B i, j P i, j+1 " B i, j-l P i,j-1 “ C i , j P i+l,j 

" E i/3 P i-l/j+l " E i+1 , j -1 P i+l,j-l " F i,j P i+1 , j+1 

- 4 D. . P* . - 4 D. . P* , . 

1/3 1+1/3 1/3 1-1/3 

n D* n t>* 


- C. , . P* . . 

1-1/3 1-1/3 


- F. , . , P* . , 

l-l, 3-1 l-l, 3-1 


— r - r ■ - I J • J f J * 

I 4 D i,j p i » j « - 4 D i.j n , a - a + D i - i , 3 n - 2 ,% 

- B i,J < + D i-1,3 * D t-t, j P i-l,j-l - 4 P U,3 

+ D,,, . P} ,43., .!>« , . .. + D. . . . ?* . , 


+ Ri+l,j ~i+2,j + D i+1 , j P i+l , j+1 + D i+1, j P i+l,j> 

+ D i/j-l P i-l,j-l 


- 4 D i+1, 




J. / J. **- 1 ^ / J * J. / J • 

. , p* . + D . , p* 

D 4-f j”-l x - j. , j 

+ D i,j-1 P i / j -2 - 4 D i, j-l P i,j-1 + D i, j+1 P £+l, j+1 
+ D i,j+1 P i-1, j+1 + D i,j+1 P i, j+2 " 4 D i,j+1 P i, j+1 


/ - 


-G. 


i/j S 


A i/j P i/j + B i , j y i, j ” B i, j-l ^i, j-l + C i, j v i, j 

C i“l / j V i~l / j + E i/j °i,j E i+1 , j-l a i+l , j-l 

4- F. . 0. . - F. , . , 0. , . . + 4 D , . L. . 

1/3 1/3 1-1/3-1 l-l, 3-1 1/3 1/3 

D i-l,j L i-l,j °i+l,j L i+l,j D i, j-l L i , j-l 


Vs " D i/ j+l L i,j+l 


Note that all terms in S and G except A . in S . . and 

, . : i / 3 f,3 

“A- • P. • in G. . involve first guess pattern information 

; x • J i / J : + / J ■ 

which is constant during the analysis . 
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° e bitten as 

S . . p * 

x ° ^ ' (G i-J + «i (j > . 0 

” H i,j' let us group tcgether 

at <^ch poi nt . 6 COeff icien ts of p, 

i#j " ^*i _7 . p* 

1,y i-2, j + (-C. . 

J i-1, i ~ 4 D. 

♦ ,-c. . 1<3 ' 4 

+ “i«, pj.. . + Di 'i " 4 « 

i# j 


J P* 

1 W,j 


1+1,j ^ +2 '0 + + D 


, 1+ l-j 

+ D i/j+1 ) P* 

+ " 4 D. . . " 

1 ' J D j_ i + i ) P * 

* ( ' F i j + D. . 

3 ?' 3 W ^, 3+1 

Di ‘ 1,j + %j-l> . 

1 1 >3-1 

- - “ 4 D i H-i) P* 

,J 1 i/j-1 


+ (-F. 

I-.,- ^o-i 1 _, , 

1 J-/ J 

+ f-B. . 

Vj*l- 4 Pi . 

l 0 


De fine; 


X 'J ' p* 

+ (-E. ' J 1 1/j-J 

1+1 + D 

+ B i i , p* , i+1 -i + D i i ,) p* 

i,3 ~ 2 + P i,j + 2 i+1 ' j 

Xi 'i B c i,i - < < D + . n f r 

1+i ' 3 

1,3 ( °i,i + D } 

J 1 '3+1 J 

+ °i+l r j + D. 

E 3 1,j+1 

i+1 0 + °i , + D. 

° - : i+lpj+1 


Y. . 

*0 

2. . 

A /j 


r\ . 

If J 


N ° te th at x, y - 

, ' 2 and R have a 

nalvs i e- a 


the analysis. 


Val ue duri 


ng 
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Then 


- H. , - D, 


. P< 


. - X. 


. P* 


. - X. • P* 


i+1 


x , 3 'i-1 , j i-l/j i x r 3 

+ D i+l,j P i+2 / j + R i-l/j P i-l,j+l Y i/j P + 


+ z. 


i,j P i+l,j+l + " "i/j- 


+ Z, 


P* 


- Y. 


P i+l,j~l i,j-l i,j-2 


D i,j+1 



,j [ 1 - 4 ] 

, j+1 

. P* . , 
1 1 , 3-1 

P* . , _ 

1 , 3+2 
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The minimization equation [1.3] is solved by simul- 
taneous over-relaxation. The matrices S. . and G. . may 

i/3 

be computed initially except for the A. . term and will 

i/3 

not change throughout the analysis. Matrix H. . must be 

1 / J 

recomputed for every iteration of the relaxation. 

The relaxation proceeds as follows: At Point (i,j) 

the terms of the minimization equation are evaluated using 
the assembled P field for P*. In general, the equation is 
not satisfied and a residual is defined as 







[1.5] 


The superscript t is an iteration counter. The value of 

P* . is to be altered so that on the next iteration, the 
i/3 

residual will be zero, provided H. . does not change. Of 

i/3 

course, H. . will change, but if the equation is fairly 
■ i / 3 

well behaved, repetition of the procedure should lead to 
convergence on the correct solution. 



P * T+1 _ 
i/j 



H i 


■ ; 



[ 1 . 6 ] 


Subtracting [1.6] from [1.5], 



i •% 



p* T +1 
1/3 


) 



[1.7] 
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and 

p* _ p* ^ 

i/j i/j 

Convergence can be hastened by increasing the correction 
term in [1.7]. The factor by which it is increased is 
called the over-relaxation coefficient. Equation [1.7] 
becomes 



p* _ p* T 

i/j i/j 


ALFA 



[ 1 . 8 ] 


ALFA must be between 1 and 2 . 

One iteration consists of making the correction 
[1.8] at every grid point. Iterations are repeated until 
the maximum residual is less than a specified convergence 
criterion. The resulting P* field is the solution of 
equation [1.3] . 
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D . Calculating The Resultant Weight 

In solving the minimization equation, we started with 

an assembled guess P. . with weight A. . . After the equa- 

1/3 i / 3 

tion is solved, the, quantity A^ ^ has been altered by adding 

a quantity, A. .P. . , which we shall call ambient information, 

i / 3 1 / 3 

The resultant value is 


A* ,P* . 
1/3 i/3 


= A. .P. • 

i/3 i/3 


A^ .P* . 
1/3 1/3 


[1.9] 


The resultant weight A* ^ will be useful in re-evaluat- 
ing the weights of the data. Holl and Mendenhall (1972) 
show that A| j could be obtained by inverting a large matrix. 
For large grids, it is not computationally feasible to do so, 
and Holl suggests a practical alternative. According to 

[1.9], a small change in the value of P, . would have the 

i/3 

following effect on the outcome of- solving the minimization 
equation: 


A* . 
i/3 


6P* , 
i/3 


= A. 


i/3 


<SP , . 

i/3 


[x- 103 


By perturbing the assembled P . .at grid points 

r - : ~. -.-I-/ 3 : , ,, 

which had data and re-solving the minimization equation 

the response in Pt . is found. From [I. 10] 
*i/3 




6P . . 
1/3 

6P* . 
1/3 
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The re-solution of [1.33 is done only in the vicinity of 

the point in question to obtain a value for 6P* . . 

1 f J 


E, Re-evaluating The Data Weights 

The quality of each observation is judged according to 
the effect that its removal has on the analysis. If an 
observation causes a change that is contrary to the guess 
field, to the differential properties of the guess field and 
to any other available observations, then its removal will 
cause a substantial change in the analysis. We have little 
confidence in such a report and would like to reduce its weight 

At the end of each cycle, the weight of each report 
is re-evaluated. For a particular report, the analysis 
value with the report removed is called the background 
analysis value, Pg, and the resultant analysis weight with 
the report removed is called the background analysis weight, A D 
To determine the impact of removing the report, A B and P B are 
computed : 



A* 


if 3 



p B = 


- A P ) /%, 
n n ' B 


[I. 11] 


The starred values are the result of the solution of [1.3] 
on the present cycle. A is the weight which the report in 
question had during the present cycle. 
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It was mentioned on page 1-2 that each weight is viewed 


as the inverse of a variance. Accordingly, the variance 

of the report is where A q is the weight the report had 

o 

originally. The variance of the background analysis value 
for this report is . Therefore, the expected difference 

a 

° 11 1/2 
between the report and its background value is + •£ ) < 

A o B 

The actual difference is P - P^. If the actual difference 
exceeds the expected difference, we shall reduce the weight 
of the report proportionately. 


Defining 


= 


< P n - P B> ' 

JL + 1 




which through manipulation and substitution of equations 
[I. 11] converts to ; 


.2 . V A l,i> < P n - P Li > 

‘ "< A i, j'VV (A i,:fV 


271 

The re-evaluated weight is defined as A „ = sr 

■’ nR l+X 2 : 

I ) ■ ' ' 

This results in: ; 

2 1 

implies actual error > expected error AA R <i 


X* >1 
2 


=1 implies actual error = expected error .VA. ^-2 
o 

X <1 implies actual error < expected error .\A nR >£ 


I -.17 


Notice that on any cycle, every data point may have its 
original weight restored, even if it had been reduced pre- 
viously. In this way, a report that causes a large change 
in the analysis may have full effect if it is supported by 
data nearby. 
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F. 


Program Description - 63 x 63 Grid Version 


1 . PCT 

The calling arguments are described in detail in 
the comments of the program. All the arrays are variably 
dimensioned, using the dimensions M and N provided in the 
calling arguments. the date-time group is provided through 
common block /DTG/ , a title for plotting, a contour starting 
point and a contour interval are in common block /INFO/ and 
sense switch variables are passed in /ISW/. 

A random-access file TAPE 9 is used for temporary 
storage and must be declared on the PROGRAM card of the calling 
program. The writing of some arrays on the random-access 
file for later retrieval allows their use as work arrays. 
Immediately after opening the file, the data list, the I 
and J data location lists, the initial data weight list and 
the initial weight field of the first guess are all written 
on the random-access file so that these arrays can be used 
in subroutine BKGRND. Since the same array is frequently used 
to hold two different fields, two names, separated by a space, 
make up the names of these arrays. Of course, the space is 
ignored by Fortran. The two names are interpreted as one, 
but this convention helps in reading the listing. 

After BKGRND computes matrices S and G (see page 1-10) , 
they are written on TAPE9 and the data-related arrays are 
read back in. The weights of the Laplacian field (D) are 
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also written on TAPE9 so that the array can be used later 
as a work array. 

DO loop 100 is the main loop controlling the number 
of cycles to be made through the program. A cycle consists of 
assembly (Section I- B), solving Equation [1.3] (Section I-C) , 
computing the resultant weight field (Section I-D) and re- 
evaluating the data weights (Section I-E) . Before doing the 
assembly, the original weights of the first-guess field are 
read from TAPE 9 . During assembly, this field is altered by 
adding the data weights. The latest solution field of P, the 
field being analyzed, is also read from TAPE9 and used as the 
guess field for this assembly. Except for the first cycle, 
the root-mean- square difference between the observations that 
were accepted on that cycle and the resulting analysis of 
that cycle is computed. 

After calling AS SMB L on the last analysis cycle only, 
subroutine PLTDAT (see Appendix) writes the data list on the 
plot file. The data list is rewritten on TAPE9 because, in 
ASSMBL, gross errors were flagged by setting the last bit 
of the data word. The last bit of ali good data is cleared. 
The assembled weight field A is written on TAPE9 as a record 
named OLDA, so that A can be used as a work array. This 
field is called "OLDA" only to distinguish it from the re- 
sultant weight which will be computed in REVALWT . Matrices 
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G and S (see page 1-10) are read from TAPE9 and subroutine 
UPDATE adds A to S and A*P to G. Next, the current guess 
field P is written on TAPE 9 and named OLDP. Weight field 
D is read from TAPE9 and subroutine BLEND solves Equation 
[1.3], resulting in a new analysis field which is written on 
TAPE9 and called NEWP. On the last analysis pass only, 
the analyzed field is passed through an optional short wave 
filter and subroutine PLOT (see Appendix) writes the analyzed 
.field on the plot file and PRT makes a printer map of it. 
Depending upon external sense switch settings the analysis 
field may also be written to the disk using the Fleet Numerical 
Weather Central (FNWC) random access routine ZRANDIO (see 
Appendix) . 

The current data weight is packed into the left half of 
each word in DWT and the original data weight in the right 
half. The original weight is needed for the weight re-eval- 
uation, Next, the guess field that was used in the call to 
BLEND on this cycle (OLDP) is packed with the latest analysis 
field (NEWP). Array DATA is used to hold these fields and to 
pass them to subroutine REVALWT. The assembled weight field 
OLD A is read from TAPE 9 and passed to REVALWT . After REVALWT 
computes new Values of DWT, arrays DATA, AI and AJ are read 
from TAPES for the next cycle. At the end of PCT, file TAPE 9 
is closed. 


2 . 


BKGRND 


Matrices S, G, X, Y, Z and R (see pages I- 10 and I- 11) 
are computed and returned. There is no problem in interpreting 
the code with the aid of the comments. 

3. ASSMBL 

The data is moved to the nearest grid point and the 
weighted average of the data and the guess value computed. 

Array SUM is first initialized to A*P, where A is the weight 
of the guess field and P is the guess field. Then, the data 
list is scanned and the guess field interpolated to each re- 
port location. If the interpolated value differs from the 
report by more than GROS, the report is rejected by setting the 
last bit of the report word. Otherwise, this bit is cleared. 

If the report is accepted, the difference between the inter- 
polated guess and the. report is added to the guess value at 
the nearest grid point and the result stored in array DATAM. 

The I and J grid point coordinates to which each report was 
moved are packed in the lower nineteen bits of each word of 
DATAM with the last bit corresponding to the last bit of DATA. 

The product of the data weight times the data value, 
after being moved to the nearest grid point, is added to SUM (I , J) 
and the data weight is added to A(I, J) . After all data have 
been scanned, SUM is divided by A to obtain the weighted 
average. 
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4. 


UPDATE 


The current assembled weight field A is added to 
matrix S (see page 1-10 and the BKGRND listing) and A*P is 
added to matrix G, where P is the current guess field. 

5. BLEND 

The minimization equation [1.3] is solved by simul- 
taneous over-relaxation. The method is described in detail 
on pages 1-13 and 1-14. No further discussion is warranted. 

6. revalwt \ ' 

The resultant weight at grid points which have data 
is computed as described in Section I-D. The weight of each 
datum is re-evaluated as discussed in Section I-E. A single 
pass through array DATAM is made. If the datum was not 
rejected in ASSMBL on the current cycle, the I,J coordinates 
of the point to which it was moved are unpacked from the 
lower nineteen bits of the word. . 

The re-solution of Equation [1.3] is done on a circle 
of three-grid-interval radius surrounding the point. The J 
limits of this area are put into JB (bottom of the area) and 
JT (top of the area) . The I limits for each row (seven rows 
total) go into the seven-element arrays IL (left of area) 
and IR (right of area) . The solution can be done only up 
to the third interior row. 


A perturbation of five percent of the assembled value 
P is added to P and the resulting change to matrix G is com- 
puted. The data weight before re-evaluation is saved to be 
printed out later. The re-solution over the limited area is 
the same as in BLEND. Once convergence is obtained, G is 
restored to its correct value and the resultant weight com- 
puted from the equation at the bottom of page 1-15. Then, 
the new data weight is calculated as specified on page 1-17 
and the datum, its old weight and its new weight are printed 
out. Finally, the guess field P is restored to the values 
it had over the limited area before the re-solution. 


G. Program Description - 187 x 187 Grid Version 
1. Design Philosophy 

A major consideration in the program design is the 
grxd size. Program specifications call for a 187 x 187 grid 
which means 34,969 computer words for each required array. 
Even if some data packing scheme is implemented, the result- 
ing arrays are too large and too numerous for the computer 
central memory capacity. Therefore, the program has been 
designed to process the 187 x 187 grid by partitions. Six- 
teen partitions are defined as follows: 

I Direction J Direction Partition 


3 

- 

48 

3 

- 

48 

1 

3 

~ 

48 : 

49 

- 

9 4 

2 

3 

- 

48 

95 

— 

140 

3 

3 


48 

141 

- 

186 

4 

49 

— 

94 

3 

- 

48 

5 

49 

- 

94 

49 


94 

6 

49 

- 

94 

95 

- 

140 

7 

49 


94 

141 

— - 

186 

8 

95 

- 

140 

3 


48 

9 

95 

- 

140 

49 

i — 

94 

10 

95 

- 

140 

95 

i 

140 ; 


95 

- 

140 

141 

- 

186 

12 

141 

- 

186 

3 

- 

4 8 

13 

141 


186 

49 

- 

94 

14 

141 


186 

95 

- 

140 

f is 

141 

— 

186 

141 


186 

16 
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The border points 1, 2 and 187 are used to process neigh- 
boring points, but are not, themselves, processed. To 
assure continuity along the borders of partitions, sixteen 
50 x 50 arrays are defined as follows: 

Partition array I Direction J Direction 


1 

1 

— 

50 

1 

- 

50 

2 

1 

- 

50 

47 

- 

96 

3 

1 

- 

50 

93 

- 

142 

4 • 

1 

- 

50 

139 

- 

188 

5 

47 

- 

96 

1 

- 

50 

6 

47 

- 

96 

47 

- 

96 

7 

47 

- 

96 

93 

- 

142 

8 

47 

- 

96 

139 

- 

188 

9 

93 

- 

142 

1 

- 

50 

10 

93 

- 

142 

47 

- 

96 

11 : 

93 

- 

142 

93 

- 

142 

12 

93 

- 

142 

139 

- 

188 

13 

139 

- 

188 

1 

- 

50 

14 

139 

- 

188 

47 

- 

96 

15 

139 

- 

188 

93 

- 

142 

16 

139 


188 

139 

— 

188 


The! points 188 are a repeat of points 187. This design allows 
a 46 x 46 partition to be processed by itself without loss of 
continuity. Note that each 50 x 50 array contains a 46 x 46 
partition. The 50 x 50 array will be referred to as a par- 


tition in this document. The term "array" is usually used to 
refer to 187 x 187 : arrays'. ' 



For an example of processing by partition, consider 
the program SSTHEM. This program uses nine 187 x 187 arrays 
having the names A, B, C, D, E, F, P, G and S. Central 
memory has storage providing for one partition at a time of 
each of these nine arrays. The remaining partitions must 
be stored outside of central memory. The program has been 
designed to store up to sixteen partitions of nine arrays in 
either extended core storage or on a mass storage file. A 
partition is brought into central memory when it is needed. 
During the processing of the partition the logic may require 
that the partition values be permanently altered. In this 
case, the partition is returned to storage with its new 
values in it. The program logic requires that old values of 
the array P be saved along with newly computed values . The 
old valued array is stored by partition on a separate mass 
storage file. - 

r The re-evaluation of data weights requires that a 
block of points from each of the arrays A, B, C, D, E, F, P, 
G, S and old P be in core for each data point to be processed 
A block must contain points plus five and minus five from the 
data point in both the I and H directions. The 16 partitions 
of an array are not adequate for this requirement. To 
accomplish the re-evaluation of data weights, an additional 
33 partitions have been defined. Each additional partition 
is a composite of the first 16 partitions. Just prior to 
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executing the PCT routine (described in the next numbered 
paragraph) , the subroutine ASSIGN is called to have a par- 
tition number assigned to each data report. Partition 
numbers are in the range 1 through 9 and are stored in the 
table KPART which has an entry for each data report. For 
convenience, the relative coordinates of the data point in 
the partition are determined and also stored in the KPART 
entry. For greater efficiency a table of the partitions 
that actually get assigned is built. This table is called 
I WHAT and contains an entry for each of the 49 partitions. 
The entry contains the number of reports that the partition 
is assigned to. A zero entry means that the partition is 
not assigned to any reports. 

The subroutine REVALWT is called on each pass from the 
PCT routine. REVALWT computes data weights, a partition at 
a time. Any partition with a zero entry in IWHAT is not 
processed. For each partition processed, REVALWT computes 
a new data weight for each of the reports having the par- 
tition assigned to it. REVALWT calls the subroutine EXTRCT 
to have array data brought into core. EXTRCT retrieves 
array data by partition number. For partitions 1 through 
16, EXTRCT retrieves a single 50 x 50 block of each of the 

y(. y 

required arrays. For higher partition numbers, EXTRCT re- 

y.l . ' . y y-:,y : .yy- y,'y.y-: 

trieves the appropriate 50 x 50 blocks and builds the re- 
quired composites. A composite set consists of a 50 x 50 
block of each of the required arrays . 
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2. PCT 


A subroutine interaction chart (Figure 1-2) for 
the program SSTHEM show the relationship of the PCT routine 
to a program package. 

When the mass storage option is selected, a random- 
access file is required for partition storage and the number 
of this file is the variable JF in the common block /MEM/. 

The value of JF must be set by the calling program. The 
TAPE9 file is opened by PCT, but the JF file is opened on 
a first call to the subroutine INDAT (see appendix) and, 
only when the mass storage option has been selected. 

The PCT logic flow consists, first, of calling the 
subroutine BKGRND to have the matrices S and G computed 
(see page 1-10) , the second, of executing one or more cycles 
of calls to ASSMBL, UPDATE, BLEND and REVALWT. DO loop 
100 is the main loop controlling the number of cycles to be 
made through the program. 

By calling INFIELD the first guess field and all 
necessary constant fields are generated, partitioned and stored 
away for subsequent processing. 

To implement partition processing the subroutine 
BKGRND is called sixteen times, once for each partition. 

Prior to each call, the subroutine INDAT is called to bring 
in a partition's worth of data for each of the arrays re- 
quired by BKGRND. After each call to BKGRND a call to 
OUTDAT (see appendix) stores a partition's worth of data 
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for each array computed or updated by BKGRND. In each 
cycle, the subroutines hSSMBL and UPDATE are called sixteen 
times using the routines INDAT and OUTDAT in the same 
manner. The BLEND routine is called sixteen times per 
iteration and the routines INDAT and OUTDAT, again, are 
used in the same manner. Iterations are executed until 
convergence is obtained, or until 150 iterations have 
been executed without convergence. When convergence is not 
obtained, PCT terminates the program run. In each cycle 
it is necessary to call REVALWT just once. REVALWT deter- 
mines for itself what partitions are required and retrieves 
these through the subroutine EXTRCT (see appendix) . 

3 . BKGRND 

This routine is identical to the 63 x 63 version. 

4. AS SMB L 

The current partition number is passed to ASSMBL 
in the calling sequence, and only reports having coordinates 
within the partition are processed. By the time ASSMBL has 
been called for all partitions, all reports will have been 
processed.. 

5. UPDATE 

This routine is identical to the 63 x 63 version. 
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6 . BLEND 


The minimization equation [1-3] is solved by simul- 
taneous over-relaxation for each partition with the maximum 
residual being passed back to PCT in common block /BLSTUF/. 

7 . REVALWT 

To implement partition processing, REVALWT makes 
use of the tables KPART and IWHAT built by the routine 
ASSIGN (see appendix) and located in the common block /REVAL/ 
In a main processing loop 49 partitions are candidates for 
processing. When a partition's IWHAT entry Is non-zero, all 
reports having the partition assigned are given an updated 
data weight on a single pass through array DATAM. If the 
datum was not rejected in ASSMBL on the current cycle, the 
I, J coordinates of the point to which it was moved are 
unpacked from the lower nineteen bits of the word and may 
be used for printout purposes. The report coordinates rela- 
tive to its assigned partition are unpacked form the left- 
most 40 bits of the report's KPART entry. (the right-most 
20 bits of the entry contain the assigned partition number.) 

The re-solution of Equation [1-3] is done in the 
same manner as in the 63 x 63 version. 
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8. Common Storage Areas - Peculiar to the 187 x 187 

Grid Version 

Blank Common 

Storage blocks in blank common have been arranged to 
accomodate a partition’s worth of data from each of nine 
187 x 187 data arrays in addition to certain other working 
arrays. For example, blank common used by the SSTI-IEM pro- 
gram is defined as follows: 

COMMON A ( 5 0 , 5 0 ) , B(50,50), C(50,50), D(50,50), 

E (50 , 50) , F (50 , 50) , P(50,50), DATA(3969), DI(3969), 

DJ (39 69 ) , DWT (3969 ) , DATAM(3969), OLDP (5-0 , 50) , 

SAVP (50,50) , SAVG (50,50) 

The arrays DI and DJ also hold data for the S and G matrices , 
respectively. This arrangement may not be disturbed without 
corresponding logic changes being made in the various routines 
which reference the blank common area. 

MEM Common 

The labeled common MEM contains cells and arrays 
defined as follows: 

Name Description 

MEMTYP This; cell is Set from a data card field. 

0 means use the extended core storage 
option. 1 means use the mass storage 
option. 
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This cell contains the number of partitions 
into which 187 x 187 arrays are divided. 

It is set to the value 16. 

This cell contains the file number of the 
storage file used in mass storage mode. 

This cell contains flags indicating the 
arrays to be brought in from storage. The 
right-most nine octal digits are used. 
Arrays 1 through 9 are assigned in order 
from left to right. An octal "1" means 
bring in a partition of the array; an octal 
"0" means don’t. 

This cell contains flags indicating 
the arrays to be written to storage. 

The flags are assigned similar to those 
of ICODE. 

This is a scratch array large enough to 
hold a partition's worth of data from 
any 187 x 187 array. 


OPTION Common 


The labeled common OPTION contains cells defined 
as follows: 


Name 


INTTYP 


KKODE 


IVEC* 


Description 


This cell is set from a data card field 
to indicate that input field data are in 
a 63 x 63 grid or a 187 x 187 grid. "0" 
means 187 x 187; "1" means 63 x 63. 

I. i • ; .. .. 

." ’ ;• ■ i : " '.'.7; •' .. . .. 

This cell contains flags indicating the 
partitions of the field arrays to be 
displayed by the routine PRT. The right- 
most sixteen octal digits of the word are 
used. Partitions 1 through sixteen are 
assigned in order from right to left. An 
octal "1" means display? an octal "0" means 
don't. (Set from card input . ) : 

This cell is set to "0" by programs 
working with scalar fields, and to "1" 
by programs working with vector fields. 


NPART 

JF 

ICODE 

JCODE 
SC (2500) 
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BLSTUF Common 


The labeled common BLSTUF is used by PCT to pass 
parameters to BLEND as follows: 


Name Description 


MM2 and NM2 These cells contain variables used by 
BLEND to limit the number of points 
processed in a partition. 

ALPHA This is the alpha variable used by 

by BLEND in the over-relaxation process 

RMAX This cell contains the current value of 

the maximum residual obtained in the 
BLEND convergence process. 


REVAL Common 


The labeled commxnon REVAL contains the following 

arrays: 


Name Description 

KPART (500) This array contains an entry for each 
data report. Each entry contains a 
partition number and a set of coordinates 
for the report. Starting from the left , 
the first 20 bits contain the coordinate 
.'I, the second 20 bits contain the coor- 
dinate J and the last 20 bits contain 
the partition number of the partition 
assigned to the report . The coordinates 
are relative to the partition areh. 

IWHAT This array contains an entry for each of 

49 partitions . Each entry contains the 
number of reports to which the partition 
is assigned, y 
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SECTION II. VECTOR WIND ANALYSIS USING THE PATTERN- 
CONSERVING TECHNIQUE 

A. Introduction 


The pattern-conserving technique described in Section I 
is used to analyze a scalar variable. In this section, we 
will concentrate on those aspects which are peculiar to the 
wind problem. 

The most essential feature of the pattern-conserving 
technique is that, while fitting new data, it tends to retain 
certain differential properties of the first-guess field. 

For scalar analysis, we were only concerned with gradients 
and the Laplacian. The wind, being a vector, complicates the 
problem slightly. Some of the properties we would like to 
conserve; e.g., vorticity and divergence, involve both scalar 
components. We must analyze both components simultaneously. 

The differential properties that we choose to conserve 
are the gradients of each wind component in eight directions 
from each grid point, the vorticity and the divergence. The 
same method is used here as in the scalar analysis, the main 
difference being that two minimization equations rather than 
one must be solved simultaneously. 


The equations are considerably 


simplified by using the 


staggered grid illustrated by Figure II-l and defining the 


divergence, vorticity and gradients as in Table II-l and 
Figure I 1-2. ■ This arrangement causes certain matrices to be 


tridiagonal . 
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The re-evaluation of the data weights has been simpli- 
fied for the wind analysis because of the increased complex- 
ity of the program and because wind reports have more in- 
herent variance than most scalar variables. 

B. Assembly 

The assembly of data to grid points is done in the same 
way as in the scalar analysis. The wind components are moved 
to the nearest grid point along the gradient of the guess 
field, and a weighted average is computed at each grid point 
for each component. Since the grid is staggered (Figure 
II-l) , the components of an observation may be moved to dif - 
ferent grid points ; i.e. , to grid points with different indices 

For the wind analysis, we need a field of assembly weights 
(A) for each component (AU and AV) . Both AU and AV are 
initialized as the class weight of the first guess. When 
a U component is assembled to a grid point (I, J) , the data 
weight of the report is added to AU (I , J) . The same data 
weight applies to both components of a report. Similarly, 
the V component is assembled to a grid point (I , J) and the 
data weight is added to A V ( I , J ) . 
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c. 


Minimizing the Deviations 


In the scalar analysis, we wanted to conserve the 
gradients and the Laplacian. With the regular grid, the 
finite difference expressions for the gradients and Laplacian 
did not provide good horizontal coupling among the grid points, 
so we included in the integral to be minimized the gradients 
and Laplacian at surrounding ' grid points. In the case of the 
wind analysis, the more complex differential properties and the 
staggered grid extend the influence of the computations at a 
grid point further than in the scalar analysis, and it is not 
necessary to add the contributions at the surrounding points. 
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TABLE II-l : PCT VECTOR CONSTRAINTS 


Constraint 

m = Variable being analyzed (assembled value) 


Weight 

A, 

l,m 


m = Variable being analyzed (assembled value) 
d l,jn = divergence = 3u/, x + 3v/ 3y 


1 ,m 


= u, , , m - u, m + v n , . - v n 
l+i .m 1 ,m 1 ,m+ 1 l^m 


D 


l,m 


l l,m = vor ticity = «lv/ 3x r.3u/ ?y 


^l.m 


v l,m v l-i,m u l,m + u l,m-i 
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We shall minimize the following integral: 



I * // [ 


(a) 

~ u l,m- 


M 

+ A l,m (v Lm ' v l,n’ 2 


(c) 

+ D 1 ,m ^ u l+i ,m u l,m + 

vt — vf - 

1 , xn + 1 1 , m 

(d) 

+ Q l,m (v l,m “ v l- i ,m ' 

ut + u? - 

l,m l,m-i 

(e) 

+ E l,m {u l-i ,n+.r u l,m 

• e l,m )2 

(f) 

/V 

+ E l,m (v i-.,n+. ' v ! ,m 

- 

(g) 

+ F l,m (u !,ir.+ : " u ! ,m " 

f ) 2 
l,m ; 

(h) 

/s 

+ F (v* — v? — 

A 2 
f l,m> 

(i) 

+ G l,m (u l+i ,m+i " u *,m 

. 2 

?! ,m 

<J) 

' ■ l •' ‘ ' . ' ' • ■ ; 

: ] ■' •; : : V - ■' . 

+ G l,r/ V l i i,m+i v l,m 

" , 2 

(k) 

; + - u X,m - 

h l n ) 2 - 

JL $ m 

(D 

+ " v l,m * 



] dxdy 


[II. 1] 
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The superscript (*) indicates the values we seek. The 
differential properties of the first guess field are defined 


in Table 11-1/ and a depiction of the u component minimization 
stencil as it relates to the u terms of equation [II. 1] Is 
given in Figure II-2 . 

To minimize the integral we take the first variation with 

respect to uf and with respect to v* , yielding the 
1 ,m 1 f m 

following two equations: 


61 

6u *l,m 

= //t 

A 1 ™ (u* _ - 
1 ,m 1 ,m 

' u l,m } 
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- 
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l,m 

* 

1+1# m 

v i 

» m~ h l,: 

) ] dxdy s ! fc 
m - 


— . . " » .X ■■ . ■■ ■ ■ T~ l — 11 ... .il 

( ) = constraint from equation [II. 1] 

*-■: = ■ U grid points- 
x = V grid points 

o = difference _ 

V 1 = divergence : v 
= vorticity 
■ — = gradient 

FIGURE II-2 : U COMPONENT MINIMIZATION STENCIL 
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In equation [II. 2] group terms involving 1) u* m ; 2) u* at 
surrounding points; 3) v* and 4) everything else. 


► ^ 

1 ,m 


[II. 4] 


J/ ((A l,m + D l,ra + Q l,ra + E l,m + F l,m + G l,rc + H l,m> 


l,m ' 


' + ('• D. - H, ) U* . 4 (-Q, ) U* 

l,m l,m l+i ,m ' v l,m 1 ,m-i 

• + (-E, „) u* + 


F, nj 

1 / m 1 

, + 
, m + 1 

(-G, ) 

l,m 

1+ i , m+ 

+ (-0, j 

1, m 

v* , 

+ Q, 

V? 

l,m+i 

1, m 

l-i, m 


J l, m 


A l.n °l,m + D l,m d l,m + ° 1 ,„, <»l im + E l,m e l , m + F l,m f l, 
■ + g l,m + H l,m !; J .n-J dxd y " 0 


m 


Group [II. 3] similarly: 


> n 

l,m 


In. 5] 


A A A /N A 

+ D l,m + 0 i,m + E l,m + F l,n + G l,m + H l,m>' V ‘l,m 


l,m 


L + (-G, )V* , + (- H, ) V* 

I. m l+i, m+i l,m l+i, m 


Y l,m {+ (D l,m ” Q l,m ) u l,m~ D l,m u l+i,m + Q l,m u l,in-i 
^ f A l* m V l,m +D l,in d l,m ®l,m q l,m + E 1 ,m e l,m + E 1 ,m ^l,m 


’l,m 


L+ G i m m + H, h, ] dxdy = 0 
L i,m J l,m l,m l,m 


Note that all terms in S and Z except A. in S, and 

c l,m l,m 

- A. u, in Z , involve first guess information which is 
1 ,m 1 ,m 1 ,n 

constant during the analysis. Similar conditions hold for 


S and Z. 
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Equations [II. 4] and [II. 5] can be written in matrix form: 


Si m U* + X, + Y, + Z, = 0 
' — ’ 1 / m ■ — ’1 f in —“X / m ■ 1 / m 


[II. 6] 


S n v* + X n + Y.. + Z, = 0 

=l,m — ==1 , m =l,m =l,m 


[II. 7] 


These equations must be solved simultaneously. The method of 
solution used is Liebmann over-relaxation. Using a first- 
guess for u* and v* , equation [II. 6] is, in general, not 
satisfied. A residual is defined bv: 


S. m u* L + X. m + Y. m + Z, m = R 

■ — ’1 f m ■■ — — 1 f m f m ■^-X f m 


[II.8] 


The superscript t is the iteration counter. We wish to find a 
next guess at u* such that the residual is zero, if the values 
at surrounding points do not change. 


S n u* T + 1 + X n + Y, + Z, = 0 
=l,m— =1 ,m =l,m =l,m 


[II. 9] 


Subtracting [II. 9] from [II. 8], 


I 1>m (H* T - H* T+1 ) = * 


, xT + 1 +T R 

and u* = u* 


Si 

=1 ,m 


[11-10] 
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Convergence is more rapid if the correction in [11-10] is 
exaggerated: 

u* T+1 = u* T - ALFA . [11-11] 

=1 ,m 

The over-relaxation coefficient ALFA must be between 1 and 2. 

At a particular grid point, u* is corrected by equation 
[II. 11] and v* is then corrected in an analogous way. In 
computing R from equation [II. 8] or from the analogous equa- 
tion in v*, the latest estimate of both u* and v* at surround- 
ing points is used. Some of them have been changed on the 
current iteration and some have not. 

During each iteration through the grid, the maximum 
residual is checked. When it becomes less than a prescribed 
convergence criterion, the equations are considered solved. 
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D. 


Re-evaluating the Data Weights 


In the scalar analysis, the resultant weight of the 
analysis was found by perturbing each guess value and re- 
solving the minimization equation in the neighborhood of the 
grid point. The resultant weight was used to re-evaluate 
the data weights. This procedure presents a considerable 
computational burden. It accounts for the majority of the 
computation in the scalar analysis. 

In the case of the wind analysis, such a procedure would 
be considerably more • time-comsuming since it would have to 
be done for the u and v components. Also, agreement with 
individual wind observations is less crucial than for most 
scalar observations, since wind reports are so sensitive to 
local effects, ship motions, elevation effects and many other 
problems. For these reasons, it is not considered worthwhile 
to solve for the resultant weights. 

The validity of wind reports is judged according to the 
vector difference between the reported wind and the analyzed 
wind. The analyzed wind is obtained by interpolation from 
the analysis fields. The root-mean-square difference is 
computed and averaged, for all the observations that were accepted 
on the current scan as a diagnostic output. If the report 
differs in vector magnitude from the analysis by more than 
the expected difference, its weight is re-evaluated. The 
expected difference is defined as the square root of the class 
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variance assigned to the report initially, which is the inverse 
of the original data weight. Define: 



A 


o 

where W is the nth report, W & is the interpolated 

analyzed wind, and A q is the original report weight. 

2 

If X is greater than 1, which implies actual error is 
greater than expected error, the report weight is computed as: 

2A 

^ = Q 

1 + X 2 

2 

If X is less than 1, the report is assigned the weight A q 
even if its weight was previously reduced. 
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E. Program Description - 63 x 63 Grid Version 
1. PCTWND 

All the arrays have variable dimensions. Common block 
/DTG/ provides the date-time group, /INFO/ the ident informa- 
tion and /ISW/ the switch settings. Random-access file TAPE 9 
must be declared on the PROGRAM card of the calling program 
and is used for temporary storage space within PCTWND. After 
TAPE9 is opened, the data location lists, the data weight 
list and the initial weight of the first-guess field are 
written on it. These arrays are used in the call to sub- 
routine BKGRND, which computes matrices S, ZU and Z V. In 
the notation of Section II-C, a ( A ) referred to the v wind 

/N 

component. Since S is the same as S, no distinction needs to 

/\ 

be made. ZU is used in the program for Z and ZV for Z. 

The convention of assigning two names to an array and 
separating them by a space is used here as it was in PCT. 

Array AI S holds the I coordinate data location list ini- 
tially, but returns matrix S from BKGRND. The matrices listed 
above are written on TAPE9 and their arrays are refilled with 
their original fields. Weight fields E and F are stored on 
TAPE 9 so they can be used as work arrays later. DO loop 100 
is the main loop controlling the number of cycles to be made 
through the program. At the beginning of each cycle, the 
original field of first-guess weights is read from TAPE 9 . 
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Subroutine ASSMBL adds the data weights to the guess weights, 
moves the data to grid points and computes the weighted aver- 
age at each grid point. After calling ASSMBL on the last 
analysis cycle only, subroutine PLTWIND (see Appendix) writes 
the data list on the plot file. 

Next, the matrices which BKGRND computed are read 
from TAPE 9 . Subroutine UPDATE adds A to S , -A*U to ZU and 
-A*V to ZV. Arrays E and F were used as work arrays in 
ASSMBL, so their fields are refilled from TAPE9. Subroutine 
BLEND solves Equations [II. 6] and [II. 7], resulting in the 
analysis fields U and V. On the last pass only, U and V 
are written on the plot file by subroutine PLOT and a 
printer map is made by PRT. Depending on external sense 
switch settings the analysis fields may also be written to 
the disk using the FNWC random access routine ZRANDIO. 
Finally, the data weights are re-evaluated by REVALWT as 
described in Section II-D, and a vector root-mean-square dif- 
ference between the observations that were accepted on that 
cycle and the resulting analysis of that cycle is computed. 

i. '"; 2;, BKGRND 

Matrices S, ZU and ZV are computed as indicated on 
page II-9. The comments adequately describe each step. 
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3. AS SMB L 


All data is moved to the nearest grid point and the 
weighted average of all data and the guess field are computed 
for each grid point. The assembled weight at each grid point 
is the guess weight plus the weight of each datum which was 
moved to that grid point. 

Arrays SUMU and SUMV are used to add up the product 
of weights times the data for each grid point. They are 
initialized to the guess value times its weight. Then, all 
the data is scanned and moved to the nearest grid point. Since 
the grid is staggered (See Figure II-l) , the J coordinate of 
the u component report and the I coordinate of the v component 
report are decreased by .5. Then, the guess u and v are 
interpolated to the adjusted report location and if a gross 
error has not occurred, the difference between the report 
and the interpolated guess is added to the guess value at 
the nearest grid point. Because of the staggered grid, a u 
report may be moved to a grid point which is labeled differently 
from the point to which the v report is moved. 

The assembled weight field is the single array A, 
but two arrays are needed: one for u and one for v. These 

two fields are packed into array A. 

Gross errors are rejected by setting the last bit of 
the DATA word. For good data, this bit is cleared. 

Finally, the weighted average of u and v are computed 
for each -grid point. 
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4 . UPDATE 


Matrix S applies to both the u and v equations, but 
the terms A*U and A*V have been left Out. UPDATE adds them 
in and two matrices result. These two are packed into array 

S. Also, A*U is subtracted from ZU and A*V from ZV. 

5. blend 

The two minimization equations, [IX. 6] and [II. 7] , 
are solved as described in detail on page 11-10. No fur- 
ther discussion is needed. 

6. REVALWT 

Data weight re-evaluation is much simpler for wind 
analysis than for scalar analysis. The explanation in Section 
II-D is followed closely and the comments in the listing 
suffice to explain the code. 
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F . Program Description - 187 x 187 Grid Version 

1. Design Philosophy 

The vector wind analysis program has been designed 
to process the 187 x 187 grid by partitions in a fashion si- 
milar to that of the scalar case described under Section I-G-l 
However, in the wind analysis case, the 16 partitions of an 
array are adequate for the re-evaluation of data weights. As 
a result, the subroutines ASSIGN and EXTRCT (described in the 
Appendix) are not required. 

2 . PCTWND 

Except for partition processing and a different 
arrangement of storage arrays, the description of PCTWND 
under Section II-E-1, also describes the 187 x 187 grid 
version of PCTWND. The partition processing implemented 
in PCTWND is similar to that implemented in PCT described 
under Section I-G-2. 

Except for blank common, the common descriptions 
under Section I-G-8 also apply to the vector wind analysis 
case. 

For the vector case blank common has been arranged 
as follows: 
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COMMON A ( 5 0 , 5 0) , DQ(5Q,50), U(50,50), V(50,50), 
ZU (50 , 50) , ZV (50 , 50) , S(50,50), DATA(63,63), 

AI (63 , 63) , A J (63,63) , DWT(63,63) 


In the 187 x 187 grid case, the matrices E, F, G 
and H have been reduced to scalars, and these scalar values 
are set by a DATA statement in PCTWND. 

3 . BKGRND 

The most significant difference between the 187 x 187 
grid and the 63 x 63 grid cases is that in the 187 x 187 case 
BKGRND does not call the subroutine. DIVERG. DIVERG is called 
from the main routine, WINDHEM. 

4 . AS SMB L 

Except for partition processing, the description of 
ASSMBL under Section II-E-4 also describes the 187 x 187 grid 
version of ASSMBL. 

The current partition number is passed to ASSMBL 
in the calling sequence, and only those reports having coor- 
dinates within the partition are processed. By the time AS SMB 
has been called for all partitions, all reports will have been 
processed. 
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5 . UPDATE 


This routine is identical to the 63 x 63 version. 

6 . BLEND 

The two minimization equations, [IT. 6] and [II. 7], 
are solved for each partition with the maximum residual 
being passed back to PCTWND in common block /BLSTUF/ . 

7 . REVALWT 

To implement partition processing, REVALWT calls 
INDATV (see Appendix) to have a partition's worth of each 
required data array brought into central memory. This is 
done for each data report processed. 


SECTION III: MASS STRUCTURE LINEAR TRANSFORMATIONS 

A. Introduction 


If an analysis of upper-air pressure heights and tem- 
peratures is to be used in initializing a forecast model, 
it is desirable for the heights and temperatures at grid 
points to be consistent with the hydrostatic equation. It 
will be shown in Section III-B that the heights, temperatures 
and layer stabilities can be interrelated through various 
linear transforms. It turns out that to close the equation 
set it is necessary to also specify a single height parameter 
and a single temperature parameter. 

The vertical organization of height and temperature levels 
and stability layers to be used in the mass structure analysis 
is shown in Figure 4. The stability parameter used here is 
defined as: 


_ E i Sln9 

P Sp 


[III.l] 


Other definitions are possible, as discussed by Holl et ad 
(1963). This definition makes 0 linear in p (k= R/Cp) , 
which is consistent with pseudo-adiabatic diagrams. 

A limitation of this technique is that a must be assumed 
to be constant in each of the layers labeled 1 - 10 in Figure 
III-l. If a serious: departure from this condition occurs in 
a layer, the temperature above the layer will depart from 
the reported temperature , but will agree hydrostatically 
with the analyzed heights. 
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Millibars 





B. Derivation of the Transformations 


The hydrostatic equation states 


dz 

dp 


1 _ 

gp 


RT 

pg 


[III ,2] 


or 


T = - P2. dz 
R dp* 

Potential temperature is defined as 


— K 


0 = T (—-) 
o 


v.’here fc. = R/ 


Cp 


[III. 3] 


[III .4] 


Therefore , 


0 


g dz 1 -k 

-|dpp 


(p 0 ) K . 


Defining n = 1-K/ 

In 0 = In ~ + In '(-—■) + nlnp +KlnP 

is Qp O 


[III. 5] 


[III. 6] 


and 


= a [ln( _dz ]+ a 

dp dp dp p 


[I II. 7] 


and 


dz dln0 
dp dp 


d 2 z 


d P 


n dz 
2 p dp * 


[III . 8] 


Substituting for p from the hydrostatic equation 
into [III .1] 


dz 2 dln0 
0 ■?: dp p dp 


[III. 9] 
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I 


and 


dlnQ _ a 
dp dz 2 

* rr 


[III .10] 


Substituting [III .10] into [III. 8] 


d z g dz 

d P 2 p dp ‘ 


[III .11] 


Multiplying by p 


n d z , n-1 dz 

> 1 — ~ + np tt s 


a T) — 2 

— p 1 

9 r 


[III .12] 


Let us define 


n dz 

X = - P dp ' 


[III. 13] 


Then 


~ - OP 


n-i dz n d“z 


[III .14] 


From [I!- .12] and [III. 14] 


ZX = Z p'l- 2 . 

dp g p 


[III .15] 


Integrating within the layer characterized by constant stability 


X + C, =-- - i-p p n 1 + C. 
A 1 g n-1 


[III .16] 
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Defining M = 

„H-1 


X ** 


- a 


g(n-D * 

Integrating [III. 13] 
Z « C, - / dp. 


+ M. 


[III. 17] 


[III .18] 


Substituting [III..17] into [III. 18] 


Z =C 3 + / tgTH=T) P' 1 - HP"") d P 


or 


K 


z=c 3 + iTOT ln P - M r^r pl_n + c 4 


and finally 

Z == N* + M* p 
where 

N* = c 3 + C 4 

and 


g< 


In p 


M* = - 


M 

l-n ‘ 


[III. 19] 


[III. 20] 


[III. 21] 


Equation [III • 21] is the basic equation of the method 
relating pressure height to stability . We also need a 
relationship between the temperature and the stability. 
Taking the first derivative of equation [III. 21] with 
respect to pressure, 


dz ... k- 1 a 

= M* K p - — , 

dp c gKp 


[ 111 . 22 ] 


III-5 


Substituting [III. 22] into [III. 3] 


T = - E£ ( M * K p K_1 

K. 


a 


gKp 


), 


[III. 23] 


or 


T = 


a 

Rk 


- M * p 


K 


[III. 24] 


Equation [III. 24] is the basic equation of the method relating 
temperature to stability . 

Repeating equation [III. 21] 


Z - N* + M* p K - In p. 

gK 


[III. 21] 


Equations [III. 21] and [III. 24] are the two model equa- 
tions we need. They apply to each of the ten layers in Figure 
III-l. The N*, M* and a in the ten layers make a total of 
30 unknowns. Applying equation [III. 21] to each mandatory 
level gives us twelve equations: 


N, + M. P 
11 a 

N + M P 
W 1 1 b 


K 


k 


Va 

a iA> 


= z 


- Z, 


[III .26] 


N io + M io P k 
N 10 + M 10 P 1 


K 


C lA 

a io e i 


= z, 


= Z. 


where 


3 = In p and the subscript (*) has been omitted. 
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Requiring continuity of height at the interface 
levels between each layer leads to nine more equations: 

N 1 - N 2 + P m K (M 1 - M 2> ' <°1 - °2> = 0 HU. 27] 

• . 0 • 

• .# • 

• ' • ' • 

N 9 - »„ + (Mg - M 3 q) - S u (o 9 - a 9 ) * 0 

Requiring continuity of temperature at the interface 
levels gives, from equation [III. 24], the remaining nine 
equations : 

a m (M 1 - M 2 ) - (a 1 - a 2 ) = 0 [III. 28] 



The vector Z_' is composed of the twelve mandatory level heights 

and 18 zeroes. The vector C is the 30-element column vector 
N 

(M) , where the ten elements of N, and M and correspond to 
a 
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the ten layers. Equation [III. 29] is written out in Figure 
III-2 . This can be represented in a partitioned form as 


A A 


*1 X 2 X 3 


N 

■ ■ i 

Z 



/s 



H 1 H 2 H 3 


M 

'55' 

0 

G 1 G 2 G 3 


a 


0 


It can be solved using Gauss elimination, giving 




[III. 30] 


Which in partitioned form can be represented as 



N 


F 1 F 2 F 3 


z 

A 





M 

— . 

E 1 E 2 E 3 


0 

a 


D 1 D 2 D 3 


0 


[III. 31] 


In the analysis, we need a transformation to get stabi- 
lities from heights. That transformation is part of matrix g ^ , 

A A 

namely a = !Z. ■ 

We will also need a transform back to heights. For that 
problem, our set of 30 equations contains 32 unknowns (10 Ns, 

10 Ms and 12 Zs) , Two of the unknowns will have to be given 
in order to close the set. 

The obvious choice for one of them is the 1000 mb height, 


since 


more data is available at the surface than in the upper 


air. Choosing the Second parameter is more difficult. Since 
the temperature will be computed from the heights, the second 
parameter might be chosen as a reference for the temperature 
profile. We chose the thickness of the 1000 - 300 mb layer. 
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( z 

Define the 12-element column vector Z_ =1 H 

where Z & is the 1000 mb height and H is the 1000 - 300 mb 

thickness. We need the transformation = D Z, 

The last ten rows of matrix g are the first twelve 

columns of the last ten rows of g ^ (D^ of [III. 31]). The 

first two rows of g are 

100000000000 

- 100000010000 . 

-1 

Matrix g may be obtained by Gauss elimination, and the 
heights can be recovered using Z_ - g ^ 

Let us repeat the sequence of operations. At grid 
points, matrix g^ is used to convert twelve mandatory level 
heights to ten layer stabilities. The stabilities are limited 
to greater than zero and less than a maximum value if necessary 
Then matrix g ^ is used to compute the mandatory level heights 
at the grid points. 

The temperatures at the grid points can now be computed 

by simply substituting a and M* , which are submatrices of 

-1 ~ 

B , into [III. 24]. In matrix form, T = Q Z_ where T com- 

. ’ ; • . \ A 

prises the twelve mandatory level temperatures, Z the twelve 
mandatory level heights and Q the matrix of equation [III. 24], 
namely 

D 1 /Rk - p K | 1 . [III. 32] 
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The temperatures and heights at the layer interface 
levels can be obtained by simply changing the pressures 
in the matrix coefficients. By proper substitution, the 
following matrix equation can be formulated: 

F = £ Z 

where 

T. 

- 5 ( Zi } 

Vectors and Z^ are the nine temperatures and heights at 
the layer interfaces . 
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SECTION IV: THE SEA-SURFACE TEMPERATURE ANALYSIS 

A. Introduction 

The Pattern Conservation Technique (PCT) is applied 
to analyze the sea-surface temperature. For Details of 
the PCT, see Section I. 

Since reports of the sea temperature tend to be less 
reliable than most atmospheric data, they are given a rela- 
tively low initial weight. Actually, the effect of lowering 
the data weight is achieved by increasing the weights of 
the guess field shape parameters. The gradients in four 
directions from each point and the Laplacian of the guess 
field may have weights as high as 0.1, compared to the in- 
itial data weight of 1.0. The shape parameter weights are 
proportional to the gradients of the terrain field, which 
are on file TAPE3 . Reducing these weights where terrain 
gradients occur has the effect of inhibiting the spreading 
of the influence of data across land. 

Since the first guess at the sea temperature (the pre- 
vious analysis) is usually fairly accurate, and since the 
data quality is not so good, the PCT is very well suited for 
this analysis. 

A brief description of the major program elements not 
covered in Section I follows. 
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B. 


SSTHEM 


The arrays necessary to do the analysis are in blank 
common. Those arrays which pertain to data must be dimen- 
sioned at least as big as the number of grid points being 
used (3,969 for the 63 x 63 grid), because they are used 
in the PCT routine to hold gridded fields. 

Common block /DTG/ holds the date-time group in two 
forms: one for identifying the data and one for display 

purposes. Common block /INFO/ contains a title for display 
and a contour origin and contour interval for plotting. 

The guess field is normally the previous analysis, but 
if sense switch 1 is "on", the guess field is read from a 
history tape. System routine SSWTCH checks the sense switch. 
If the sense switch is "on", subroutine READNED searches the 
history tape for the desired field, which is described by 
array IPT. READNED is described in the Appendix. If switch 
3 is on, the 12 hour old guess field is read using ZRANDIO. 

If switch 3 is off, the guess field is read from the last 
analysis output tape (TAPE5) using RDRITE . If switch 4 is 
on, the output fields are written using ZRANDIO. If switch 5 
is on, the output fields are written to TAPE 8 using the sub- 
routine PLOT. If switch 6 is on, data lists, analysis sta- 
tistics and a contour map of the final analysis are produced. 

Subroutine READATA or ADPSST extracts the sea temperature 
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reports from the data files. Next, bogus cards are read 
from the input file until an end-of-file is encountered. 

Each card is one report. Columns 1-10 of each card should 
contain the bogus Centigrade temperature report in normal 
floating point format. If the decimal point is punched, 
the numbers may appear anywhere within the ten columns. If 
not, columns 9 and 10 will be interpreted as the tenths and 
hundredths digits. The same rule applies to the location 
and weight parameters. Columns 11-30 should contain the grid 
point location of the bogus report. The I coordinate goes 
in Columns 11-20; the J coordinate in Columns 21-30. Columns 
31-40 should hold the weight which the user wishes the bogus 
data to have initially. 

Next, the terrain gradients in four directions from each 
point are read from TAPE3 and statement function WT computes 
the weights of the corresponding gradients of the guess field. 
The weight of the Laplacian is the average of the four grad- 
ient weights. The weight of the guess value is .01. 

The number of passes to be made through the PCT program 
is specified by the variable NOP AS . The convergence criterion 
to be used in the PCT program is CONV. If a report differs 


from the guess value interpolated to the report location 


by more than GROS, the report is summarily rejected for 


that pass. All reports are rechecked at the beginning of 
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each pass. 

All that remains is to call the PCT subroutine. See 
Section I and the comments in the program for the details 
of the calling arguments. 

C. READATA/ADPSST 

The routine READATA searches the file containing all 
types of reports (referred to as the "spot data") for ship 
reports of sea-surface temperature. The temperatures are 
put into array DATA; and X and Y grid point coordinates of 
the report location go into arrays DI and DJ, respectively; 
the initial data weights go into array DWT; and the total 
number of reports found is put into the variable NOREP. 

The routine ADPSST 'searches the FNWC disk resident 
ADPFILE for ship records containing sea surface temperature 
observations. Knowing the FNWC data format the observations 
are decoded and put into the arrays, described above. 

The formats for both methods of accessing data are 
peculiar to the FNWC data base and further elaboration and 
explanation are not deemed necessary. 
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SECTION V. TEMPERATURE AND SEA-LEVEL PRESSURE GUESS 
A . Introduction 

Two different methods can be used to generate the tem- 
perature and sea- level pressure guess fields. The method 
used is dependent upon whether an analysis/forecast cycle 
is available or whether one is producing the first analysis 
in a cycle. If no primitive equation forecasts are avail- 
able, the first-guess fields for analyzing the sea- level 
pressure and upper-air temperature fields are computed by 
advecting the previous analyses with half of the smoothed 
500 mb geostrophic wind. Free-slip boundary conditions are 
used in the advection computation. Program MAKGS generates 
this type of first-guess field. 

It has been found by operational experience that if a 
12-hour forecast verifying at the analysis time is available, 
these produce better first guess fields than advection of 
the last analysis, especially in the upper air. This method 
is used after the initial analysis is produced in the anal- 
ysis/prediction/analysis cycle. Program GSFCST generates 
this type of first-guess fields. 
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B. 


Advective First-Guess 


MAKGS 


If sense switch 1 has been turned on prior to executing 
MAKGS , it is assumed that no previous analyses have been 
made and that file TAPE1 is the history file of the fields 
to be advected. If switch 2 is on, current date-time-group 
fields are read from the disc using ZRANDIO and no advection 
is performed. This option is used to initialize the 12-level 
analysis from the 10-levels available at FNWC . The two 
new levels (950 and 900 mb) are obtained from a linear 
interpolation in ln(p) using the 1000 and 850 mb levels. 

If switch 3 is on the 12-hour old input fields are read 
using ZRANDIO. If switch 3 is off the 12-hour old fields 
are read from the last analysis output tape (TAPE5 ) using 
RDRITE . If switch 4 is on the output fields are written 
to the disc using ZRANDIO. If switch 5 is on the output 
fields are written to TAPE6 using the subroutine PLOT. If 
switch 6 is on the output fields ar'e contour mapped. 

First, the 12-hour-old 500 mb height field is read 
and smoothed very heavily by subroutine FILTR. (See Sec- 
tion III-4-C of the report on the Primitive Equation Fore- 
cast Model.) The geostrophic wind, is then computed from 
the smoothed field. 

The 12-hour-old sea-level pressure field is then read. 
Subroutine ADVEGT advects the old pressure field forward 
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for twelve hours, using one-hour time steps. The resulting 
guess field for the sea- level pressure analysis is smoothed 
lightly and output in a manner determined by the sense 
switches. 

Finally, the 12-hour-old upper-air temperature fields 
are read and advected for twelve hours with the same wind 
field used before. The resulting guess fields are smoothed 
and output in a similar manner. 
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C. ADVECT 

Field X is advected for the number of hours specified 
in the variable ISTEP, using the wind field supplied in 
arrays U and V. The winds are modified along the boundary 
by subroutine BCWND (see below) . The field being advected 
is first modified along the boundary to remove any gradient 
normal to the boundary. Array DX2 is filled with the dis- 
tance in meters of two grid intervals, taking proper account 
of the map factor. 

The advection is begun with a forward time step, fol- 
lowed by successive centered time steps, until the specified 
number of one-hour steps has been made. Simple centered 
space differencing is used in the advection. After each 
time step, subroutine BC is called to ensure that no gradient 
normal to the lateral boundaries can build up. 

D. BC 

The lateral boundary points of field X are made equal 
to the second interior point, removing any gradient normal 
to the boundary across the first interior grid row. The 
corner points are filled by averaging the three surrounding 
points. 


V-4 



E . BCWND 

Fields U and V are changed so there is no wind normal 
to the boundary on the first interior grid row and no 
gradient of the component tangent to the boundary in the 
direction perpendicular to the boundary. 

F . Prognostic First-Guess GSFCST 

The program GSFCST extracts from the PE forecast out- 
put tape (TAPE10) the appropriate 12-hour forecast surface 
pressure and upper-air temperature fields. If sense switch 
4 is on the output fields are written to the disc using 
ZRANDIO . Is switch 5 is on the output fields are written 
to TAPE 6 using PLOT. If switch 6 is on the output fields 
are contour mapped. 
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SECTION VI. THE HEMISPHERIC SEA-LEVEL PRESSURE ANALYSIS 
A. Introduction 

The same techniques used in SSTHEM are also used in 
the sea level pressure analysis (PSI-IEM) . The guess field 
is prepared by program MAKGS or GSFCST prior to the running 
of PSHEM. 

If switch 3 is on, this guess field and data are read 
using ZRANDIO; if off, the field is read from TAPE6 using 
RDRITE . If switch 1 is on, the input data is read using 
READATA from a "spot data" TAPE1. Similar to SSTHEM, if 
switch 4 is on the output field is written using ZRANDIO; 
with switch 5 on the output is written to TAPE8 using PLOT. 
If switch 6 is on, data - sample and contour maps are printed. 

The pattern conservation technique is used to preserve 
the shape of the guess field according to weight fields A 
through F. The data is fitted moreclosely in SPHEM than in 
SSTHEM by lowering the weights of the shape parameters from 
a maximum of .1 to .01. The terrain gradients in file TAPE3 
are used to decrease the weights of the corresponding shape 
parameters where the terrain is uneven, just as in SSTHEM. 

In nature, mountain ridges act as obstacles to the migration 
of surface pressure patterns. An observation on one side 
of. such a mountain ridge should not materially affect the 
analysis at grid points on the opposite side . The reduction 
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of the weights simulates this effect. 

Bogus data is read from data cards in the same way as 
in SSTHEM. 

B. PSHEM 

All the necessary arrays are dimensioned in blank 
common. The arrays for holding information about the data 
are dimensioned larger than those holding gridded fields 
because of the large number of surface pressure observa- 
tions that are available. 

Common block /DTG/ holds the date-time group for iden- 
tification of data (IDT) and for display ( JDT) . Common 
block /INFO/ contains a title and a contour origin and 
interval for plotting. Common block /ISW/ contains the 
sense switch indicators. 

The input field and data are read in a manner depending 
upon how the sense switches are set. Then the bogus reports 
are read from data cards until no more cards remain. Each 
card contains one report. Columns 1-10 hold the reported 
pressure in mb; 11-20 hold the I coordinate of the report 
location; 21-30 the J coordinate; 31-40 have the weight to 
be assigned to the report. All four parameters should be 
in F10. 2 format. If no decimal point is punched, it is 
assumed to be just left of the second column from the right 
side of each ten-column field. A punched decimal point 
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overrides this rule. 

Statement function WT computes the weights of the 
shape parameters based on the gradients of the terrain field 
in four directions from each point, which are read from 
file TAPE3. 

The PCT program is called to perform the analysis. 

The method of output is determined by the sense switch 
settings . 

C. RE ADATA/ADP P S 

Similar to the description in SSTHEM, READATA searches 
a "spot data" tape for sea level pressure reports from ship 
and land stations. ADPPS extracts the same type of reports 
from the appropriate FNWC disc resident operational ADPFILE 
records. Whichever way the data is accessed, the report 
value is put into the array DATA, the locations into DI and 
DJ and the weight into DWT. 

The data is assigned a weight value of 1.0 for land 
reports and 0.7 for ship reports. There is a problem in 
decoding pressure reports, since the hundreds digit is left 
off. A report more than 50 cannot be assumed to be less 
than 1,000 mb. Pressure greater than 1,050 mb does occur; 

1 . f ■ : ' - ■ V • . fv . 

likewise, a repqrt of 49 could be 949 mb. To resolve this 
problem, the guess pressure field, which is passed to the 
data extraction routine as an argument, is interpolated to 


the report location whenever the report is between 40 and 
65. If the guess is less than 1,000 mb, then the report 
is assumed to be less than 1,000 mb, also. 

Notice that only reports at the exact synoptic time 
are accepted because the pressure field can change quite 
rapidly . 
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SECTION VII. UPPER-AIR GEOPOTENTIAL HEIGHT AND 
TEMPERATURE ANALYSIS 

A. Introduction 

The mass structure analysis procedure will be described 
in four main parts: the virtual temperature analysis, the 

computation of the height guess fields, the height-tempera- 
ture analysis and statospheric field extrapolation. 

The virtual temperature analysis uses the PCT and guess 
fields at each pressure level, prepared by MAKGS or GSFCST 
(See Section V) . Then the sea-level pressure analysis 
(already done by PSHEM) is hydrostatically converted to a 
1,000 mb height field using the 1,000 mb virtual temperature 
analysis. The analyzed virtual temperatures are then used 
to compute the other heights . 

These height fields are used as the guess fields in 
the height analysis program. The PCT is again used, but 
this time, in addition to using the height data from radio- 
sondes, all available wind data (aircraft reports, pibals, 
radiosonde reports) and satellite wind reports are used to 
compute geostrophic height gradients. The gradient data 
replaces the gradients of the guess field at those points 
where wind data is available. These gradients are given mu cn 
higher weights than those computed from the guess field. 

After the height analysis is made, matrix D (see page 
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III-10) is multiplied by the twelve-element column vector 
Z (the twelve mandatory level heights) to obtain the sta- 
bilities in each of the ten layers. (See Figure 4.) If 
the heights imply a hydrostatically unstable layer, the 
stability of that layer is made slightly positive and ma- 
trix (see Page III-10) is used to compute the corrected 

heights. Finally, matrix Q (see page ITI-10) yields the 
final virtual temperature analysis. 

This technique can be described as a mass structure 
technique in that the final temperatures and heights are 
hydrostatically consistent and stable. 

There are at least three advantages of this method 
over the analysis of stabilities as done by FNWC: 

1. Incomplete soundings can be used 

2. Wind data is used directly 

3. Satellite-observed temperature data can be 
used without a geopotential reference 

A potential problem area in this analysis is the 
matter of a geopotential reference. At present, sea level 
is the basic reference level. The 1,000 mb height is com- 
puted hydrostatically from the sea-level pressure and the 
1,000 mb temperature, assuming that the layer between sea 
level and 1,000 mb has a mean temperature equal to that of 
the 1,000 mb level. Then the 1,000 mb height is used as 
the reference level for all the other heights. The advantage 
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is that the many sea- level pressure reports are utilized. 

The disadvantage is that many of these reports are from 
high-altitude stations and the sea-level pressures have 
little meaning. If a great deal of satellite temperature 
data becomes available, it is worthwhile to consider using 
an upper level, where satellite-observed temperatures are 
most reliable, as the reference level. 

Finally, 50, 30 and 10 mb height and temperature fields 
are produced using an empirical extrapolation method which 
builds upon the 100 mb analyzed height and temperature 
fields . 
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B . Temperature Analysis - TMPANAL 


The PCT is used to analyze radiosonde virtual tem- 
perature observations at the twelve mandatory pressure 
levels. The guess field has been prepared by program 
MAKGS or GSFCST. (See Section V) 

All the necessary arrays are in blank common and are 
passed as arguments to the PCT program. Common block 
/DTG/ holds the date-time group for data identification 
(IDT) and for display (JDT) . Common block /INFO/ contains 
a display- coded title and a contour origin and contour 
interval for plotting. Common block /ISW/ contains sense 
switch setting information. 

If sense switch 1 is on, the input data is read from 
a "spot data" tape. If switch 3 is on, the first guess 
fields and data are read using ZRANDIO. If switch 3 is off, 
the input fields are read from TAPE 6 using RDRITE . If 
switch 4 is on, the output fields are written using ZRANDIO; 
with switch 5 on, the output is written to TAPE6 using PLOT. 
If switch 6 is on, data samples and contour maps are printed 
The weights of the shape parameters (gradients and 
Laplacian) are set to a constant value so the terrain does 
not affect the analysis. 

For each of the twelve mandatory levels, the first 
guess field and data are read in a manner specified by the 
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switch settings and PCT does the analysis. For each level, 
the INPUT file is checked for bogus data. Bogus reports 
should be ordered by pressure level, one report per card, 
with the levels separated by blank cards. Any levels for 
which no bogus data is provided should be represented by a 
blank card. 

The bogus data should occupy the first ten columns of 
each card in F10.2 format. The I and J coordinate location 
and initial data weight go into columns 11 through 40 in 
F10. 2 format. 

The temperature analyses are output in a manner spec- 
ified by the switch settings. 

C. READ AT A/ AD P UAT 

The data file is searched for radiosonde temperature 
observations for date-time IDT and level LVL. The source 
of the input data is determined by switch settings. 

If a dewpoint temperature is reported, the virtual 
temperature is computed by subroutine VIRT. If not, the 
temperature report is accepted as the virtual temperature. 
The entire analysis is done in degrees Centigrade. 


D. 


Height First-Guess - HTGUES 


The program HTGUES calculates the first-guess height 
field. If sense switch 3 is on, the sea-level pressure and 
1,000 mb temperature fields are read from the disc using 
ZRANDIO. If switch 3 is off, the input fields are read 
from TAPE8 using RDRITE. The height of the 1,000 mb surface 
above sea level is computed by integrating the hydrostatic 
equation, assuming that the 1,000 mb temperature represents 
the average temperature for the layer between sea level and 
1,000 mb. If switch 4 is on, the resulting 1,000 mb height 
field is written to the disc using ZRANDIO. If switch 5 is 
on, the fields are output to TAPE7 using PLOT. Then the 
other mandatory level temperatures are read and the other 
heights computed similarly. If switch 6 is on, the resulting 
output height fields are contour' mapped. 
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E . lie ight-Temp erature Analysis -x'liT^NAL 

The heights of the mandatory pressure levels are ana- 
lyzed with the PCT program, using guess fields computed by 
HTGUES and radiosonde and satellite retrieved height reports. 
Wind data from radiosonde reports, pibals , aircraft and 
satellite derived reports are used to compute height gradients 
which are treated as data in the PCT program. 

After the analysis, the stabilities of the ten layers 
shown in Figure 4 of Section III are computed using matrix 
D. {See page III-10.) Any instabilities are removed and 
the virtual temperatures at the twelve mandatory levels are 
computed using matrix Q. {See page III-10 . ) 

All the necessary arrays are in blank common and are 
passed as arguments to the subroutines. Common block /DTG/ 
holds the date-time group for data identification (IDT) and 
for display {JDT) . Common block /INFO/ has a title, a con- 
tour origin and a contour interval for display. Common 
block /LVL/ contains the number (1-12) of the mandatory 
level being analyzed and /ISW/ contains the sense switch 
indicators. 

If. sense switch 1 is on, the data is read from a "spot 
data" tape (TAPE1) . If switch 3 is on, the input fields 
and data are read using 2RANDI0. If switch 3 is off, the 
input fields are read from TAPE7 using PLOT. Switch 4 on 
prompts the output fields to be written to the disc using 
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Z RAND 10. If switch 5 is on, the output height fields are 
written to TAPE8 and the temperature fields to TAPE4 using 
PLOT. If switch 6 is on, data samples, analysis statistics 
and contour maps are produced. 

For each of the twelve levels, the guess field and the 
height data are read in a manner specified by the sense 
switch settings. Then the input file is checked for bogus 
data in the same manner as described in Section IV-B. 

Bogus reports should be ordered by pressure level, one 
report per card, with the levels separated by blank cards. 
Any levels for which no bogus data is provided must be rep- 
resented by a blank card. The bogus height should occupy 
the first ten columns of each card in F10.2 format. The I 
and J coordinate location and initial data weight go into 
columns 11 through 40 in F10.2 format. 

The PCT subroutine does the analysis and writes the 
output in a manner specified by the switch settings. 
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F. READATA/ADPUAH 


These subroutines are exactly like those described in 
Section VII-C, except that heights are read instead of 
temperatures. Reports are given an initial data weight 
of .0025. 

G. RDWNDS 

Subroutine BKGRND (see Section I) , which is part of 
the PCT, has been slightly altered to call RDWNDS after 
computing the gradients of the guess field. RDWNDS calls 
subroutines AIREP, PIBAL, RASONDE and ADPSAT, which read 
wind data from aircraft, pilot balloon radiosonde and sat- 
ellite wind reports, respectively, and compute corresponding 
height gradients. The gradients of the guess field are in 
arrays XMU (positive J direction) and XNU (positive I direc- 
tion) when RDWNDS is called. Subroutines AIREP, PIBAL, 
RASONDE and ADPSAT put the height gradients into the corre- 
sponding arrays YNU and YMU at the points corresponding to 
the lower left corners of the grid squares within which the 
wind reports lie. A look at the gradient definitions in 
Table 1-1 on page 1-5 will clarify this point. If more than 
one report is available for a grid point, the one implying a 
height gradient closest to the gradient of the guess field 
is accepted. The weights on the gradient terms (E and C) 
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are increased at points where wind data is used. To deter- 
mine whether data has been accepted at a grid point, the 
weight value B at the point is compared to B{1,1). If it 
is different, then data was used. This assumes that data will 
not be available at point (1,1), a safe assumption, since 
data south of the Equator is rejected. 

H. AIREP/ADPAIRP 

Aircraft-reported wind data are read in a manner spec- 
ified by the switch settings. Reports within six hours of 
the time of the analysis are accepted, but their initial 
weights are reduced in proportion to their age. 

The flight level of the reporting aircraft is used to 
find which mandatory level the report is nearest. The 
standard heights of the pressure levels are used for this 
comparison. The report is accepted unchanged at the closest 
level. 

Once the wind direction and speed are unpacked, the u 
and v components on the polar stereographic grid are com- 
puted by subroutine DDFF2UV. The X and Y gradients are 
computed geostrophically as: 

XGRAD = £ • V* AX 

g 

YGRAD ” - • U* AY 

g 
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Here, f is the Coriolis parameter? g the acceleration of 
gravity; AX and AY the ' grid distances in the X and Y di- 
rections. If the gradients differ from those of the guess 
field by more than a specified amount, the report is 
rejected. 

The value of the weight of the gradient in the Y 
direction is checked against the value at point (1,1). 

If they are not equal, then data must have already been 
utilized at the point in question. In this case, the 
gradients computed from the previous report (already in 
YNU and YMU) are compared to the gradients of the guess 
field (XNU and XMU) . The differences are called OLDIFI 
and OLDIFJ . A similar comparison is made between the grad- 
ients just computed from the report being considered and 
the guess gradients. The differences are called NEWDIFI 
and NEWDIFJ. If the sum of the new differences is less 
than that of the old ones, the new report is accepted. 

The weights on the gradient, data are computed from the 
age of the report. A report at the time of the analysis 
gets a weight of .002; a six-hour-old report has a weight 
of .0015. The guess gradients were assigned a weight of 

.ooi. f ^v; 

Array IPLT was passed as an argument from RDWNDS . 

The u and v components of the wind report and its I and J 
coordinates are packed into a single word of IPLT . 
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Winds from pilot balloon observations are read and 
treated as those from aircraft reports were treated in 
AIREP with the source being specified by a switch setting. 
The only significant differences between AIREP and PIBAL, 
except those resulting from the slightly different data 
format, are that only reports for the exact time of the 
analysis are accepted and the gradient weights are all .002. 

J . RASONDE/ADPUAW 

Winds from radiosonde observations are read and treated 
exactly as in PIBAL, The only differences are due to the 
data formats. 

K. ADPSAT 

Winds derived from satellite cloud vectors are read. 

The reports include a pressure level estimate and are 
included in the analysis at the closest analysis level to 
this estimate. The computed gradients are given a value 
of .0015. 

L. PLTWND ■ ■■!• 

This subroutine is called by RDWNDS with array IPLT 
as an argument. Array IPLT contains the components and 
the I, J location of all the wind reports which were used 
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in RDWNDS. A twenty-word identification package is added 
as detailed in the description of subroutine PLOT in the 
Appendix and an unformatted record containing IPLT and the 
identification array is written on the plot file. The fifth 
word of the identification array is set to 7 HWNDATA to 
identify tne record as packed wind data. 

M. TEMP 12 

After the PCT program finishes the analysis on all 
twelve pressure levels, TEMP12 is called to check for hydro- 
static instability and to compute the final temperature 
analysis from the heights. As discussed on page III-10 , 
the 12 x 12 matrix D, when multiplied by the twelve-element 
column vector Cthe twelve mandatory level heights at a 
point) yields the stabilities of the ten layers shown in 
Figure 3-1. If a negative stability is computed, a message 
is printed and the stability is changed to a slightly pos- 
itive value. Excessively high stabilities are reduced. 

Then matrix DINV (the inverse of D) is multiplied by the 
twelve element vector consisting of the 1,000 mb height, 
the 1,000 to 300 mb thickness and the ten stabilities to 
give the corrected heights. See page III-10 for details. 

Finally, matrix Q is multiplied by the corrected heights 
to produce the final temperature analyses. The explanation 
of this transformation is on page III-10 . 
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N. 


Stratospheric Height-Temperature 
Extrapolation - STRATQ 


The method of obtaining stratospheric height and tem- 
perature fields is by vertical extrapolation as presented 
by Lea (1961) . The extrapolated values are given by equa- 
tions of the form 


Z~A o +A i z ievel-l +A 2 T level-l 


and 


T=A 0 +A.,Z, - i+ArT n . , . 

3 4 level-1 5 level-1 


In this manner the 50 mb fields are extrapolated from the 
100 mb height and temperature, 30 mb from the 50 mb and 
10 mb from the 30 mb. The coefficients which are a function 
of latitude and month were obtained from empirical studies 
of selected rawinsonde stations. The empirical constants 
are defined in the block data routine STRTCON. After the 
fields have been produced they are filtered once by a long 
wave filter. 
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SECTION VIII. UPPER-AIR WIND ANALYSIS 


A. Introduction 

The application of the Pattern Conservation Technique 
(PCT) to wind analysis is discussed in Section II. The 
discussion in this section is limited to the peculiarities 
of the upper-air wind analysis program. The technique is 
used to analyze the winds at the twelve standard pressure 
levels . . 

Recall that the • PCT wind analysis fits the wind data, 
while conserving the gradients in four directions, the 
vorticity and the divergence of the guess field. Each of 
these differential properties has an associated field of 
weights. In the upper-air wind analysis, the surface 
terrain is not used to compute these weights. 

Since the analyzed wind fields are used to initialize 
the Primitive Equation Forecast Model (PE HEM) , special 
attention should be given to controlling the divergence of 
the final winds. Excessive initial divergence, which is 
not reflected in the initial temperature and geopotential 
fields, causes large- amplitude, high-frequency oscillations 
to develop early in the forecast. This phenomenon is usually 
referred to as geostrophic adjustment and can seriously 
degrade the usefulness of the forecast for twelve or more 
hours, until the mass and motion fields become balanced. 
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Because the divergence is a very sensitive differential 
property of the wind field, unrealistically large values 
will usually result from any analysis scheme unless the 
divergence is explicitly controlled. The PCT is an effec- 
tive and versatile method of controlling the divergence. 

A common philosophy on this .subject is that zero is 
the best available estimate of the initial divergence. 

This reasoning led to the use of the balance equation to 
derive a wind field containing no divergence from the 
geopotential field, . While zero is certainly a better 
estimate than the large values which arise from a scalar 
analysis of each wind component, the twelve-hour forecast 
divergence from a previous run of the PEHEM is likely to 
have some skill. 

The PCT can be used to apply either philosophy. A 
slight modification is made to the basic PCT wind program 
to read the wind field whose divergence is to be conserved 
from a file, instead of computing the divergence from the 
guess field. If no divergence is desired in the analyzed 
wind, switch 1 is left off. If the twelve-hour forecast 
divergence is desired, switch 1 is set and the PE forecast 
winds are read from TAPE10 using RDRITE . 

The first-guess wind field at each of the twelve 
standard pressure levels is the geostrophic wind computed 
from the analyzed geopotential field for that level . 


OS' POOR QUALITY 
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B. 


WINDHEM 


The PCT technique is used to analyze the winds at each 
of the twelve pressure levels. If sense switch 2 is on, 
the input data is obtained from a "spot data" tape,. TAPE1. 

If switch 3 is on, the height analysis and input data are 
read from the disc using ZRANDIQ. If switch 3 is off, the 
input fields are read from TAPES using PLOT. If switch 4 
is on, the output fields are written using ZRANDIO; if 
switch 5 is on, they are written to TAPE11 using PLOT. If 
switch 6 is on, partial data lists, analysis statistics and 
contour maps are produced. 

Subro-utine WIND computes the geostrophic wind to be 
used as the first guess. Subroutines AIREP , PIBAL, RASONDE 
and ADPSAT read wind data from aircraft reports, pilot 
balloons, radiosondes and satellite wind reports at the 
proper level indicated by variable LVL, which is in a common 
block. 

Bogus data is read from the INPUT file until a zero 
data weight is found. Bogus reports for different levels 
should be separated by a blank card. When the blank card 
is read, the data weight will be zero and the reading is 
stopped. Any levels for which no bogus reports are to be 
inserted must be represented by a blank card. If no bogus 
reports are desired, twelve blank cards should be in the 
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INPUT file. 


One card is required for each bogus report. The 
direction from which the wind blows relative to north goes 
in the first ten columns, the speed in meters per second 
in the next ten columns, the I and J coordinates of the 
report position in columns 21 through 40 and the initial 
data weight in columns 41 through 50. Each ten-column 
field is read with an F10.2 format. 

Subroutine DDFF2UV converts the bogus wind data to 
u and v components. The components are packed into the 
single array DATA with u in the left half of each word and 
v in the right. The entire data list is written on the 
plot file by subroutine PLTWIND . 

Finally, the PCTWND subroutine does the analysis. 

C. WIND 

The geostrophic wind components are computed from the 
geopotential height analysis. In the computation of the 
grid distance and the Coriolis parameter, the latitude is 
not allowed to be less than twenty-five degrees. The geo- 
strophic wind law is inadequate in tropical latitudes and 
has a singularity at the Equator. 
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D. 


AIREP/ADPAIRP , PIBAL/ADPPIB , 
RASONDE/ADPUAW , ADPSAT 


. The first three subroutines are identical to those by 
the same names described in Sections VII-H, X and J , 
except that the computation of the geopotential gradients 
has been deleted and the u and v wind components are packed 
into array DATA. The locations and initial data weights 
are put into arrays AI, AJ and DWT. 

ADPSAT reads cloud vector derived satellite wind 
reports if switch 2 is on. These reports contain a pres- 
sure level estimate in centibars at which the observation 
is believed valid, an age, and a direction and speed. The 
observation is included in the data list of the pressure 
level nearest the observation. The observation is given a 
weight of 0.5 in comparison to 0.1 for a radiosonde. 

*/ 
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E. 


BKGRND 


In the PCTWND program (see Section II) subroutine 
BKGRND has been changed to call DIVERG, This subroutine 
reads the twelve-hour forecast wind field and computes its 
divergence. 1% BKGRND, array ZU is used temporarily to 
hold the computed divergence. Normally, both the diver- 
gence and vorticity fields have the weight DQ. We now 
wish to give the divergence a higher weight. To accom- 
plish this, the divergence weight is packed into the left 
half of each word in PQ and the vorticity weight into the 
right half. A minor change was necessary in subroutine 
BLEND of PCTWND to allow the divergence and vorticity to 
have different weights. 
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F. DIVERG 

A wind field for level LVL (in common block /LVL/ ) 
is read from TAPE 10 if switch 1 is on. The output for the 
initial time is skipped. The u component field is read 
and the contribution of the gradient of u to the divergence 
is computed. 

Then the v component is read and the contribution of 
the v gradient is added to the divergence field. The final 
divergence field is smoothed heavily. \> 

The divergence field is assigned the same weight as 
wind reports have. This weight field is packed with the 
vorticity weight field in array DQ. The divergence weights 
are in the left half of each word and the vorticity weights 
in the right. 

If switch 1 is off, the divergence value is set to 
zero. ■ 
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SECTION IX. SUPPORT PROGRAMS 


A. CKRAOB 

Radiosonde data and satellite retrieved radiosonde 
profiles are known to occasionally contain errors due to 
garbling in transmission or errors in the actual workup 
of the sounding. As mentioned earlier , the upper air 
analysis is constrained to hydrostatic equilibrium. The 
hydrostatic equation can be used as a check of whether 
reports at successive levels in the sounding are reason- 
able. Due to radiosonde malfunction, levels may be missing 
in the sounding and the hydrostatic equation can be used to 
approximate a reasonable value for the missing levels. The 
mass structure analysis is required to analyze heights and 
temperatures at the 950 and 900 mb levels. These are not 
mandatory levels for radiosonde reports so the observation 
must be interpolated from the profile that is available 
and written in a format the analysis programs can read. 

If sense switch 1 is off, the input data is read using 
Z RAND 10; if on, the data is read from a "spot data” tape. 
First, significant level observations are read and dupli- 
cates removed. Next mandatory level observations are read 
and duplicates are removed. Next the mandatory and signif- 
icant levels of observations reporting both are merged into 
an array in which pressure monotonically decreases (height 


increases) . Tropopause and maximum wind data are also 
inserted at the appropriate level if reported. 

The next steps are checking for hydrostatic consis- 
tency, assigning heights to significant level reports, and 
approximating missing levels if possible. 

First the temperature lapse rate is checked using 
3 degrees C/100 meters as a gross check. If the lapse rate 
exceeds this, the level is flagged as missing and the check 
continues until the top of the sounding is reached. Next 
an attempt is made to find the station level report. This 
is supposed to be the first level reported in a significant 
level report. Knowing the station elevation, one can cal- 
culate the standard station pressure. If the two values 
differ by greater than 60 mb, this is probably not the sta- 
tion level report. If not, the first reported mandatory 
ievel is used for the base height for the hydrostatic workup 
Other miscellaneous checks are performed to try to further 
verify the station level report. The checks are adequately 
described by comment cards so one can understand this logic. 
Given a surface pressure, height and temperature, one can 
calculate, using the hydrostatic equation, the next pressure 
level height knowing the pressure and temperature at that 
level. These heights are included in the significant level 
reports and changed in mandatory levels if the report does 
not appear to be consistent. Below 250 mb, the average 
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height change required to make the sounding hydrostatically 
consistent is on the order of 10 meters in a sample size of 
approximately 450 soundings. Above. 250 mb the change runs 
from about 20 meters to 50 meters. Next, missing tempera- 
tures and dewpoint levels are interpolated in ln(p) if 
possible. Wind reports from both mandatory and significant 
levels are merged with tropopause and maximum wind reports. 
Missing levels are interpolated in ln(p) if possible. 

Finally, the levels analyzed by the mass structure 
analysis are extracted and written to the disc in a format 
compatible with the analysis programs input. The reports 
are extracted as reported, if possible, or interpolated in 
ln(p) if no reports are available at the needed level. 

This applies to heights, temperatures, dewpoint depression, 
wind direction and wind speed. 


B. 


INITIAL 


The purpose of the program INITIAL is to generate the 
input file for the primitive equation model initialization. 
The routine BREAD is used to read from the appropriate tape 
files the needed initialization fields in the correct order. 
The order is as follows: 


Fields 

Description 

Units 

1-13 

Heights (1000 - 50 mb) 

meters 

14-26 

Temperatures (1000 - 50 rob) 

degrees Celsius 

27 

Land-sea- ice table 

integer 

28 

Sea level pressure 

millibars 

29 

Terrain height 

meters 

30 

Sea surface temperature 

degrees Celsius 

31 

Albedo 

percentage 

32-55 

u/v wind field pairs 
(1000 - 100 mb) 

meters/second 


The fields are written to the output file using BWRITE . 

'The heights and temperatures are those produced by 
HTANAL. The land-sea-ice table is an integer field with 1 
denoting land, 2 meaning sea and 3 meaning ice. This field 
was obtained from the FNWC operational system and is repre- 
sentative of the day-time-group of the test data. The sea 
level pressure was analyzed by PS HEM. The terrain height 
was obtained by analyzing one degree terrain data obtained 
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from Scripps Institute of Oceanography onto the grid and 
limiting the normalized terrain gradient to less than 2000 
meters per standard mesh length. This terrain is similar 
to that used by the FNWC operational primitive equation 
model. The sea surface temperature field was analyzed by 
SSTHEM. The albedo field was obtained by interpolating 
climatological fields produced by Posey and Clapp to the 
polar stereographic grid. These fields are also similar 
to those used by the FNWC operational model. The wind 
fields were analyzed by WINDHEM. 


c. 


IDENT 


The program IDENT is used to produce inventories 
of any of the output files generated by the analysis pro- 
grams that have been written by BUFFER OUT. This includes 
analyzed fields and data records. The twenty word idents 
which are part of these field are decoded and pertinent 
information printed. 

D . PLTDWA 

This program is executed to produce output that can 
be processed for display on the Versatec plotter. The 
input file is TAPE2, which can be any of the output tapes 
produced by the analysis programs. By decoding the ident 
of the program determines which type of field it has read 
and processes it for the plotter in the appropriate manner. 

E. CREATE 

An option has been included in the analysis programs to 
allow access to data and initial fields independent of the 
FNWC operating system, specifically ZRANDIO. This option is 
exercised by setting IREAD in common block/MS/ equal to one 
(1) vice zero (0) and having the sense switches set as if to 


read the data with Z RAND 10 . 

The program CREATE is executed at the beginning of 
the program execution sequence (see Section X-B) . The 
input tape (TAPE1) contains in BUFFER OUT format one tape 
per date-time-group. The tape contains an index record 
which specifies which data records and fields follow, then 
these records. An index record specifying the names of 
the data index records is next, followed by the indices. 
The index records contain the names of the individual 
records of data, one index for each specific type of data. 
The execution of CREATE transfers the input tape to disc 
in a READMS random access format on TAPE16 . When read by 
the analysis programs, the data is transferred to central 
memory in exactly the same format as a ZRANDIO read. 


SECTION X: PRODUCTION ORGANIZATION 

A. Introduction 

The execution of the analyses is performed in a 
sequence of programs executed within a single computer job. 
Each successive program writes its output to the disk with 
ZRANDIO or to a TAPE file with BUFFER OUT so the fields are 
available to succeeding programs. When the sequence is com- 
plete, one file has been produced which contains all the 
analysis output fields and another file has been generated 
in the format needed to initialize the prediction model 
(see IX-B) . 

All the programs that produce upper air fields have 
an option for producing a 10 or 12 level output. The ten 
level model does not output the 950 and 900 mb levels. 

The number of levels is specified by the number 10 or 12 
in parenthesis on the execution card. Coding is included 
to gain access to this number and use it as an indicator 
within the program. 

B. Program Execution Sequence 

The actual sequence of programs that are executed depend 
upon whether it is the first analysis execution in the day- 
time-group series or whether output from previous analyses 
or forecasts are available. 
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In the first analysis in the series the current FNWC 
analyses are used as the first guess for the normal analysis 
levels. Since 950 and 900 mb analyses are not produced at 
FNWC, these levels are interpolated in ln(p) from the 
available 1000 and 850 mb levels. Since the current analysis 
is available, it is used as the first guess instead of an in- 
ferior twelve hour advected field. 

Therefore the program sequence is as follows: 


Step 

Program 

Output 

0. 

CREATE 

Random access data and initial 
fields (optional) . 

1 . 

CKRAOB 

Checked radiosonde reports. 

2. 

SSTHEM 

SST anal. 

3. 

MKGUES 

Produce surface pressure and upper air 
temperature first guess fields from cur- 
rent FNWC fields, interpolate 950 and 
900 mb fields. 

4. 

PSHEM 

Analyze surface pressure. 

5. 

TMPANAL 

Analyze 1000 to 100 mb temperature fields 

6. 

HTGUES 

Generate upper air height analysis first 
guess from temperature anals (1000 to 
100 mb) . : 

7. 

HTANAL 

Analyze height fields and retrieve con- 
sistent temperature fields (1000 - 100 mb 
Extrapolate 50, 30 and 10 mb strato 
heights and temperatures. 

8. 

WINDHEM 

Analyze upper air wind fields (1000 - 
100 mb) . 

9. 

INITIAL 

Generate forecast initialization file. 
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Prior to the execution of the sequence the required 
fields and data must be made available. Upon completion, 
the output files are copied to tape for later use. 

In subsequent analysis cycles, the output from the pre- 
vious analysis or prognostic fields from a prediction ini- 
tialized from the previous analysis are available. In these 
studies twelve hour prognoses are used as first guess fields 
for surface pressure and upper air temperature analyses. In 
this second type of execution sequence Step 3 of the above 
sequence is replaced, with: 

3 . GSFCST Read prediction model output tape and 

write surface pressure and upper air 
temperature first guess fields. 

This implies that for each analysis sequence after the 
initial day-time-group a twelve hour forecast must be calcu- 
lated to use as first guess for the next analysis sequence. 
The first guess for the sea surface temperature analysis is 
the twelve hour old analysis field. Therefore the output 
from the last analysis cycle must also be available. 

The TAPE numbers and content involved in the analysis 
sequence are given below: 
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Tape 

Content 

Program 

TAP El 

NEDN "spot data" (not used) 
or data/field tape (optional) 

All anals read 
CREATE read 
(optional) 

TAPE 2 

NEDN field tape (not used) 

SSTHEM read 

TAPE3 

Terrain gradients 

SSTHEM read 
PSHEM read 

TAPE 4 

Retrieved temperature anals 
and strato temperature 
extrapolations (temporary) 

HTANAL write 

TAPES 

-j . 

Last analysis sequence output 

SSTHEM read 
MKGUES read 

TAPE 6 

Surface pressure and upper 
air temperature first guess 
(temporary) 

MKGUES write 
GSFCST write 
PSHEM read 
TMPANAL read 

TAPE7 

Height anal first guess 
(temporary) 

HTGUES write 
HTANAL read 

TAPE 8 

Analysis sequence output tape 

All write 
INITIAL read 

TAPE 9 

WRITMS/READMS Analysis 
scratch file 

All anals read 
and write 

TAPE 10 

12 hour old forecast output 
tape 

GSFCST read 
WINDHEM read 

TAPE 11 

Wind anal output (temporary) 

WINDHEM write 

TAPE 12 

Land/sea/ice field 

INITIAL read 

TAPE 13 

Terrain height field ' - 

INITIAL read 

TAPE 14 

Albedo field ; - it,;-' 

INITIAL read 

TAPE 15 

Prediction model initializa- 
tion fields 

INITIAL write 

TAPE 16 W;-V.- 

READMS data/initial fields 
■ ■■;■■■ (optional ) 

All anals read 
(optional) 

TAPE 17 

187x187 model READMS /WRITMS 
partition scratch file 

All 187x187 anals 
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As soon as HTANAL and WINDHEM have finished executing, 

TAPE4 and TAPEll are rewound and copied to TAPES. This 
results in the following fields and data records being 
written to TAPE 8 , the analysis sequence output tape, for 
the twelve level analysis output. 

Number Record Total Num- 

Type Field Levels Types/Level ber Records 

1 3 data/1 field 4 

1 3 data/1 field 4 

12 3 data/1 filed ' 48 

12 4 data/1 field 60 

3 1 field 3 

12 1 field 12 

3 1 field 3 

12 3 data/2 fields 60 

Total =194 


Sea surface temp. 

Surface pressure 
Upper air temp, anal 
Height anals 

Strato height extrapolation 

Temperature retrievals 
Strato temp, extrapolation 

Upper air wind anal 





C. Production Series 

The data sets which were saved from FNWC operational 
fields and data included eight successive twelve hour date- 
time groups (DTG) . The dat'i sets included appropriate ana- 
lysis and prediction fields for verification and all possible 
types of raw data that might be used. 

The first DTG will be used to initialize the twelve 
layer analysis model* This means using the first execution 
sequence as described in Section X-B. From this analysis 
cycle a twelve hour forecast will be run. This forecast will 
be used as the first guess for the next DTG. This analysis 
should have the influence of the FNWC analysis fields removed 
and forecast divergence .will be available to input to WINDIiEM. 
The execution sequence will be the second type in Section X-B. 
This DTG will be Considered the base date-time-group for the 
prediction model. Seventy-two hour forecasts of the various 
prediction models will be run from this DTG. Analyses to use 
for verification will be produced using the second type of 
execution sequence, namely running twelve hour forecasts to 
produce 'first guess fields for the next analysis DTG. This 
will be done for the remaining six date-time-groups until an 
analysis set is produced to verify the seventy- two hour forecas 
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LATLNIJ 

The latitude/longitude location of a point is trans- 
formed to I and J grid coordinates. Latitudes and Longi- 
tudes are input in degrees. The trigonometries of the prob- 
lem are pictured in the diagrams below. Triangles ABC and 
AFE are similar, so 

r R 

d D 

The unknown value is R, a distance to be measured in grid 
lengths. Angle <p is the latitude; R q is the earth’s radius. 

r = R cos'eb 
e 

d - R e (l+sin<j>) 

Since 0 is 45°. D = RED, which is the number of grid intervals 
from the Pole to the Equator. (The. standard value of RED for 
the hemispheric 63 x 63 grid is 31.205.) Therefore, 

_ RED*cos<j> r* n 

R - i + iisy • tA - u 

R is the radial distance on the map from the Pole to the 
point in question, measured in grid lengths. The longitude 
must be used to calculate I and J. 1 

The longitude is input as a positive angle in the Eastern 
Hemisphere and negative in the Western. It is converted to 


A-l 



a "mathematical style" angle, increasing counterclockwise 
from the positive I direction by subtracting 10° from it. 
Designating this angle as FAKLON and R as FACTOR, it is 
easily seen that 

I = 32 + FACTOR • cos (FAKLON) 

J « 32 + FACTOR • sin (FAKLON) 
for the 63 x 63 grid. 


A 



RED 


•CROSS-SECTION EARTH POLAR VIEW EARTH 


FIGURE A- Is 
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The I , J grid coordinate location of a point on the 
63 x 63 polar stereographic grid is transformed to latitude 
and longitude. The latitude can be solved directly from 
Equation [A.l] in the previous paragraph. 




2 2 

,RED - R , 

' 2 2 ' 

RED + R 


The longitude is easily computed by subtracting the I and J 
coordinates of the Pole from the coordinates of the point and 
using the ATAN function. A look at Figure A-l v. ill clarify 
this point. 


PLOT 


All the analysis programs use PLOT to write results to 
a disk file. A twenty-word identification is added in front 
of the array, giving the information necessary to direct the 
plot program to plot the field. Subroutine RDPLT can be used 
to read the fields written by PLOT and to remove the twenty- 
word identification group. 

The calling arguments include the array to be written 
(DATA) ", a four-word title (ITL) , a single-word display coded 
date-time-group for labeling the plot ( IDTG) , a starting 
point for computing contours (CO) , the interval between con- 
tours (Cl), a number for computing labels on contours (CL - 
not presently used) , the forecast time in hours, an integer 
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(ITAU) , the unit number of the file to which the field is to 
be written, an integer (IUNIT) , the dimension of the field 
in the X direction, an integer (M) and the dimension of the 
field in the Y direction (N) . All this information is put 
into the twenty-word identification group as shown in Table A-l . 

PLOT always writes a single unformatted binary record 
of 3,989 words. Common block /PLOT/ holds the array that is 
used for this purpose. Subroutine RDPLT also uses this common 
block. If the array being written out is smaller than 3,969 
words, there are some unused words at the end of each record. 
Since the dimensions of the field are in words seven and eight 
of the record, the proper field is easily retrieved. 

TABLE A-l: TWENTY-WORD IDENTIFICATION GROUP 

Word Content Word Content 


1 

IDTG (Display Code) 

11 

Title (Display Code) 

2 

ITAU (Integer) 

12 

Ttile (Display Code) 

3 

CO (Real) 

13 

Title (Display Code) 

4 

Cl (Real) 

14 


5 

Blank for Fields, 
4HDATA for Data 
7HWNDATA for wind 
gradient data 

15 


6 

NOREP (Integer) 

16 


7 

M (Integer) 

17 


8 

N (Integer) 

18 


9 


19 

CL (Real) / 

10 

Title (Display Code) 

20 
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PLTDAT 


The scalar analysis programs use PLTDAT to write data 
lists to a disk file to be plotted. The same file that 
holds the analysis fields is used for data. A data list is 
distinguished from a field by the fifth word of the twenty- 
word identification group described in the previous paragraph 
on PLOT and in Table A-l: For a data list, the fifth word 

is the word DATA in display code. The first record of a 
data list is 3,989 words, the first twenty of which are the 
identification group followed by the data list. The number 
of data words in the record is put into word six of the 
identification group. The 3,989-word record is followed by 
two records whose length is that of the data list with no 
identification group. These two records contain the I and 
J locations of each datum in the list. There is a serial 
correspondence among the data list, the I location list and 
the J location list. 

The calling arguments are the array containing the I 
locations, a real array (AI) , the real array containing the 
J locations (AJ) , the data list, a real array (DAT), the 
scale factor (explained below [SCALE] ) , the length of the 
data list, an integer (N) , the unit number of the disk file 
to write on (IUNIT) , the I dimension of the field being 
analyzed, an integer used to interpret the I and J locations 
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(II) and the J dimension of the field being analyzed (JJ) . 

Before writing the data list on the file, the real 
numbers in DAT are scaled, by multiplying them by the real 
number SCALE, and converted to display code. The plot pro- 
gram is designed to plot only three digits of each datum. 

SCALE should be specified so that the digits one wishes to 
plot are in the units, tens and hundreds positions after the 
multiplication. For example, if a height of 5580. is to be 
plotted, one would wish to select the first three digits. In 
this case, SCALE should be 0.1. 

PLTWIND 

The wind analysis program uses PLTWIND to write the 
wind data list to a disk file. The arguments are the same as 
for PLTDAT , except that SCALE is left out. The description of 
P LTD AT in the previous paragraph applies also to PLTWIND, 
except the data list is not scaled or converted to display 
code. Each word of the data list (DAT) contains a u and a 
v wind component, packed as real numbers. The data list, 
with the identification group before it, is written on the 
file in that form. ; 

PLTWND (In Program HTANAL) 

• Program HTANAL also writes a wind data list on its result 
file, but in order to avoid using excessive amounts of central 
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memory, the u and v components and the I and J location 
are converted to integers and packed into a single word. 

The identification group is added and a single 3,989-word 
record is written. The fifth word of the identification 
group is the display-coded word WNDATA. 

RDPLT/RDRITE 

All the analysis programs Use RDPLT or RDRITE to read 
the results of previously run programs from disk files. In 
RDPLT the data records are skipped, the identification group 
is removed and the desired array is returned to the calling 
program. Common block/PLOT/ is used to read the 3,989-word 
records. In RDRITE no distinction is made between types of 
records. One record is read for each call of the routine and 
the identification group is removed. 

ZRANDIO 7 ■■ - 7 '-’. 

With the proper switch settings, all programs have the 
capability of reading fields and data using ZRANDIO which 
is FNWC's random disk input/output routine. ZREAD, ZWRITE 
and ADPREAD are general purpose interfaces to ZRANDIO. ZREAD 
and ADPREAD are used for reading fields and observational data 
respectively while ZWRITE is used for writing fields. Unique 
labels are assigned to each record using a formula which 
combines a three character catalog number (i.e. , A01) , the 
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the date- time-group , and the forecast increment. Since the 
development and production runs of the analysis models were 
done on FNWC 1 s computers, this option of disk access was 
used whenever possible due to its greater efficiency and 
compatibility with the operating system. It would be 
necessary to implement the READMS option as described in 
Section IX-E if computers other than FNWC ' s were used. 
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INFIELD 


This routine reads input fields, partitions them, and 
either writes them to extended core storage or on a new mass 
storage file. 

The field description is specified in the calling 
parameters. The field may be read using READNED, RDRITE or 
Z RAND 10 depending upon the sense switch settings and is stored 
in its appropriate array based upon the value of IFLD. 

If ITYPE eguals "1" instead of "0", the terrain gradient 
fields are also read, partitioned and stored in their appropriate 
locations for use in the sea surface temperature or sea level 
pressure analyses. 


INDAT 


The subroutine INDAT retrieves a partition’s worth of 
data from up to nine different 187 x 187 arrays. The 187 x 
187 arrays must have already been partitioned and written 
to either mass storage or extended core storage. 

On each call to INDAT, one of a possible sixteen parti- 
tions may be specified via a calling sequence argument. This 
partition number applies to all arrays. Arrays are specified 
through the word ICODE in labeled common /MEM/ defined under 
Section I-G-8. The right-most nine octal digits of ICODE 
are used and arrays 1 through 9 are specified in the order 
left to right. A "1" digit means retrieve a partition from 
the array and a "0" digit means don’t. For example, 

ICODE = 111111111B would cause a partition's worth of data 
from all nine arrays to be retrieved. ICODE = 100100100B 
would cause a partition’s worth of data to be retrieved from 
arrays 1,4 and 7 . 

The subroutine has two entry points, INDAT and INDATV. 
When INDAT is called, the subroutine assumes that the nine 
arrays are A, B, C, D, E, F, P, G and S in that order; and 
when INDATV is called, the subroutine assumes that there are 
seven arrays in the order A, DQ, U, V, ZU, ZV and S. INDAT 
is called by scalar programs such as SSTHEM, and INDATV is 
called by the vector program WINDHEM. For an explanation of 
the arrays named here, see the corresponding program descrip- 
tions and listings . 


When retrieving data from the arrays P, U and V, the 
peripheral points of the partition are updated from points 
in the surrounding partitions where appropriate. (Peripheral 
points being those points which overlap surrounding parti- 
tions and which would not necessarily contain the latest 
updated values.) 

INDAT has been hard coded to unpack all arrays (see 
PACK and UNPACK) when processing in extended core storage 
mode. INDATV has been hard coded to unpack only the arrays 
U, V, ZU and ZV. The arrays A, DQ and S are not packed. 

Internal to INDAT, the table IN (8,16) guides the peri- 
pheral points processing. IN contains a list of eight par- 
tition numbers associated with each partition. For example, 
Partition 10 has the list 9, 11, 6, 14, 7, 5, 15 and 13 
which are the numbers of the partitions peripheral to par- 
tition 10, as illustrated here: 


13 

14 

15 


^ii;l 0 
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Note that the list is in the order left side, right 
side, bottom side, top side, lower-right corner, lower- 
left corner, upper-right corner and upper-left corner. 

Another example is Partition 8, which has the list 7, 0,4 
12, 0, 3, 0 and 11. The zero entries indicate that no 
peripheral partitions exist on the right-side, the lower-right 
corner and the upper-right corner of Partition 8, as illus- 
trated here: 


11 

12 

0 

7 

8 

0 


4 

0 


In extended core storage mode and in the scalar case, 
INDAT computes the storage address using the formula: 

I ADR = 2700* (N-l) +17* (IP-1) 

where N is an array number (1-9) and IP is a partition number 
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In the vector case INDAT makes use of the internal tables 
JADR and JFACT (see listing of INDAT) , so that the formula 
becomes: 

I ADR = JADR (N) + JFACT (N) * (IP- 1) 

where N and IP are as previously defined. 

In mass storage mode the record number is computed 
using the formula: 

IREC = 16* (N-l) + IP 

where, again, N and IP are as defined above. Prior to calling 
INDAT for the first time, the random file JF (a variable in 
labeled common /MEM/) must have been opened. This file uses 
the numerical index option. (See the description of READMS 
in the FORTRAN manual.) 
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PACK and UNPACK 


The subroutine PACK packs an array of 2500 real formatted 
numbers into the first 1700 cells of the same array, as illus- 
trated here: 


CELL 


1599 

1600 


2499 

2500 



L1701 


— 

Rl 7 01 

• 


L2500 


R2500 

• 

• 

L1701 

R1701 





• t 
• 




L2500 

R2500 






L and R indicate the left 20 bits and the right 20 bits, 
respectively, of the 40 most significant bits of a real 
formatted word. The 20 least significant bits of words 
1 - 1600 and words 1701 - 2500 are lost. 

The subroutine UNPACK unpacks the 1700 cells into the 
2500 cells using a reverse procedure. 

C18750 

The subroutine C18750 divides a 187 x 187 array into 
the 16 partitions defined under Section I.G.l. 

When the flag IOPT in labeled common /OPTION/ is set 
to " 0 " , C18750 transposes the 187 x 187 array before par- 
titioning. 

C18750 writes partitions to mass storage or to ECS, 
depending on the setting of the variable MEMTYP in the 
labeled common /MEM/. (Labeled commons are described under 
Section I.G. 8. ) - 

On each call to C18750, just one 187 x 187 array is 
processed. The array is specified via an argument in the 
calling sequence. When the argument is set to a negative 
value, the partitions are not packed, when processing under 
ECS mode. Formulas for computing ECS addresses and record 
numbers are similar to those used by INDAT . 
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C50187 


The subroutine C50187 builds a 187 x 187 array from 
the 16 partitions of the array. 

C50187 reads partitions from mass storage or from ECS 
depending on the setting of the variable MEMTYP . 

On each call to C50187 just one 187 x 187 array is 
built. The array is specified via an argument in the cal- 
ling sequence. 

C50187 always transposes the 187 x 187 array once the 
array has been built.’ 

Formulas for computing ECS addresses and record numbers 
are similar to those used by INDAT and C18750. 

TRANSPOS 

The subroutine TRANSPOS transposes an array. It is 
called by the subroutines C18750 and C50187 to have 187 x 
187 arrays transposed. 

OUTDAT 

The subroutine OUTDAT stores a partition's worth of 
data from up to nine different 187 x 187 arrays. Partitions 
are stored either on mass storage or in extended core stor- 
age depending on the setting of the variable MEMTYP. 

OUTDAT uses the word JCODE in /MEM/ in a fashion 
similar to the way that INDAT uses ICODE . 
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On each call to OUTDAT, one of a possible sixteen 
partitions may be specified via a calling sequence argument. 
This partition number applies to all arrays. 

The subroutine has two entry points , OUTDAT and 
OUTDATV. When OUTDAT is called, the subroutine assumes 
that the nine arrays are A, B, C, D, E, F, P, G and S in 
that order; and when INDATV is called, the subroutine 
assumes that there are seven arrays in the order A, DQ, U, 

V, ZU, ZV and S. Formulas for ECS addresses and record 
numbers are similar to those used by INDAT. 

LINEAR and CUBIC 

The subroutine LINEAR performs a two-dimensional linear 
interpolation of points contained in a 63 x 63 grid. 

The subroutine CUBIC calls routines to generate a two- 
dimensional interpolation of points contained in a 63 x 63 
grid using a cubic spline technique. 

The routines LINEAR and CUBIC are used to generate 
187 x 187 arrays from 63 x 63 arrays. 

ASSIGN 

The subroutine ASSIGN assigns a partition number to 
each data report. Partition numbers assigned are in the 
range 1 - 49, where the partitions numbered 17 - 49 are com- 
posites of those numbered 1 - 16 . 













































For example, Partition 17 overlaps Partitions 1 and 2 
as follows: 



and Partition 22 overlaps Partitions 2 and 6 as follows: 


i 1 



I I 


(For actual partition limits, see routine listings) . 
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ASSIGN stores a partition number into the array KPART 
which has an entry for "each report. It also stores the I,J 
coordinates of the report point relative to the partition. 
The KPART entry format is 



where I and J are the report point coordinates and IP is 
the partition number. 

As ASSIGN builds the KPART array it maintains a sum of 
reports to which each partition is assigned and stores these 
sums in the array IWHAT , which has an entry for each of 49 
possible partitions. The arrays KPART and IWHAT are con- 
tained in the labeled common /REVAL/ and are used by the 
subroutine REVALWT to retrieve partitions for data weight 
re-evaluation. 

The subroutine ASSIGN prints the contents of the IWHAT 
array just prior to returning to the calling program. 

EXTRCT 

' The subroutine EXTRCT retrieves a partition's worth of 
data from the arrays A, B, C, D, E, F, P, G, S and OLDP to 
be used by the REVALWT routine to re-evaluate data weights . 
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The required partition is specified via a calling 
sequence argument. For partitions numbered 1-16, EXTRCT 
retrieves a single partition (per array) corresponding to 
the number specified. For partitions numbered 17-49, EXTRCT 
retrieves all partitions necessary for the construction of 
the required composites. 

For all arrays except OLD? , EXTRCT retrieves partitions 
from either mass storage or extended core storage depending on 
the setting of the variable MEMTYP . EXTRCT retrieves OLDP 
partitions from the mass storage scratch file, TAPE9 . 

The internal table IN (4, 33) is used by EXTRCT to guide 
the construction of composites. Each partition having a 
partition number in the range 17-49 has an entry in the IN 
table. An entry consists of four partition numbers. For 
‘example, the entry for Partition 25 has the numbers 3, 4, 


7 and 8, which means that Partition 25 is a composite of 
partitions 3, 4, 7 and 8 as follows: 
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Another example is Partition 34 having the numbers 0, 0, 7 
and 11, which means that Partition 34 is a composite of 
Partitions 7 and 11, as follows: 


r 

i 


! 


— 34 — — 



A final example is Partition 17 having the numbers 1, 2, 0, 
0, meaning that Partition 17 is a composite of Partitions 
1 and 2, as follows: 


r — 

I.'. 

; i 

i 

i 


17 




These three examples illustrate three composite cases. 
A Case I composite is generated from quarters of four other 
partitions; a Case II composite is generated from the hori- 
zontal halves of two other partitions; and a Case III com- 
posite is generated from the vertical halves of two other 
partitions. 

Formulas for ECS addresses and record numbers are 
similar to those used by INDAT. 
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